diff --git a/doc/docs/Python_Tutorials/Eigenmode_Source.md b/doc/docs/Python_Tutorials/Eigenmode_Source.md index c89312373..627b04490 100644 --- a/doc/docs/Python_Tutorials/Eigenmode_Source.md +++ b/doc/docs/Python_Tutorials/Eigenmode_Source.md @@ -94,7 +94,9 @@ else: Note that in `EigenModeSource` as well as `get_eigenmode_coefficients`, the `direction` property must be set to `NO_DIRECTION` for a non-zero `eig_kpoint` which specifies the waveguide axis. -Additionally, we can demonstrate the eigenmode source for a rotated waveguide. The results are shown in the two figures below for the single- and multi-mode case. There is one subtlety: for mode **A** in the multi-mode case, the `bnum` parameter is set to 3 rather than 2. This is because a non-zero rotation angle breaks the symmetry in the $y$-direction which therefore precludes the use of `EVEN_Y` in `eig_parity`. Without any parity specified for the $y$-direction, the second band corresponds to *odd* modes. This is why we must select the third band which contains even modes. +### Oblique Waveguides + +The eigenmode source can also be used to launch modes in an oblique/rotated waveguide. The results are shown in the two figures below for the single- and multi-mode case. There is one subtlety: for mode **A** in the multi-mode case, the `bnum` parameter is set to 3 rather than 2. This is because a non-zero rotation angle breaks the symmetry in the $y$-direction which therefore precludes the use of `EVEN_Y` in `eig_parity`. Without any parity specified for the $y$-direction, the second band corresponds to *odd* modes. This is why we must select the third band which contains even modes. Note that an oblique waveguide leads to a breakdown in the [PML](../Perfectly_Matched_Layer.md#breakdown-of-pml-in-inhomogeneous-media). A simple workaround for mitigating the PML reflection artifacts in this case is to double the `thickness` from 1 to 2. diff --git a/doc/docs/Python_Tutorials/Optical_Forces.md b/doc/docs/Python_Tutorials/Optical_Forces.md index d4f9dfa08..340e1a394 100644 --- a/doc/docs/Python_Tutorials/Optical_Forces.md +++ b/doc/docs/Python_Tutorials/Optical_Forces.md @@ -10,7 +10,7 @@ This tutorial demonstrates Meep's ability to compute classical forces via the [M The gradient force $F$ on each waveguide arising from the evanescent coupling of the two waveguide modes can be computed analytically: -$$F=-\frac{1}{\omega}\frac{d\omega}{ds}\Bigg\vert_\vec{k}U,$$ +$$F=-\frac{1}{\omega}\frac{d\omega}{ds}\Bigg\vert_\vec{k}U$$ where $\omega$ is the mode frequency of the coupled waveguide system, $s$ is the separation distance between the parallel waveguides, $k$ is the conserved wave vector, and $U$ is the total energy of the electromagnetic fields. By convention, negative and positive values correspond to attractive and repulsive forces, respectively. For more details, see [Optics Letters, Vol. 30, pp. 3042-4, 2005](https://www.osapublishing.org/ol/abstract.cfm?uri=ol-30-22-3042). This expression has been shown to be mathematically equivalent to the Maxwell stress tensor in [Optics Express, Vol. 17, pp. 18116-35, 2009](http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-17-20-18116). We will verify this result in this tutorial. In this particular example, only the fundamental mode with odd mirror-symmetry in $y$ shows the bidirectional force. @@ -20,7 +20,7 @@ The gradient force $F$ can be computed using two different methods: (1) using MP The simulation script is in [examples/parallel-wvgs-force.py](https://github.com/NanoComp/meep/blob/master/python/examples/parallel-wvgs-force.py). The notebook is [examples/parallel-wvgs-force.ipynb](https://nbviewer.jupyter.org/github/NanoComp/meep/blob/master/python/examples/parallel-wvgs-force.ipynb). -The main component of the Meep script is the function `parallel_waveguide` which computes the force $F$ and transmitted power $P$ given the waveguide separation distance `s` and relative phase of the waveguide modes `xodd=True/False`. The eigenmode frequency $\omega$ is computed first using `get_eigenmode` in order to specify the frequency of the `add_force` and `add_flux` monitors. An [`EigenModeSource`](../Python_User_Interface.md#eigenmodesource) with `eig_match_freq=False` is then used to launch the guided mode. Alternatively, a constant-amplitude point/area source can be used to launch the mode but this is less efficient as demonstrated in [Tutorial/Eigenmode Source/Index-Guided Modes in a Ridge Waveguide](Eigenmode_Source.md#index-guided-modes-in-a-ridge-waveguide). The waveguide has width/height of $a=1$ μm and a fixed propagation wavevector of $\pi/a$. +The main component of the Meep script is the function `parallel_waveguide` which computes the force $F$ and transmitted power $P$ given the waveguide separation distance `s` and relative phase of the waveguide modes `xodd=True/False`. The eigenmode frequency $\omega$ is computed first using `get_eigenmode` in order to specify the frequency of the `add_force` and `add_flux` monitors. (Note that when `match_frequency=False`, `get_eigenmode` ignores the input `frequency` parameter.) An [`EigenModeSource`](../Python_User_Interface.md#eigenmodesource) with `eig_match_freq=False` is then used to launch the guided mode using a pulse (the `frequency` but not the `fwidth` parameter of `GaussianSource` is ignored). Alternatively, a constant-amplitude point/area source can be used to launch the mode but this is less efficient as demonstrated in [Tutorial/Eigenmode Source/Index-Guided Modes in a Ridge Waveguide](Eigenmode_Source.md#index-guided-modes-in-a-ridge-waveguide). The waveguide has width/height of $a=1$ μm and a fixed propagation wavevector of $\pi/a$. ```py import meep as mp @@ -128,7 +128,7 @@ plt.legend(loc='upper right') plt.show() ``` -The following figure shows the $E_y$ mode profile at a waveguide separation distance of 0.1 μm. This figure was generated using the [`plot2D`](../Python_User_Interface.md#data-visualization) routine and shows the source/flux region (red hatches), force monitors (blue lines), and PMLs (green hatches) surrounding the cell. From the force spectra shown above, at this separation distance the anti-symmetric mode is repulsive whereas the symmetric mode is attractive. +The following figure shows the $E_y$ mode profiles at a waveguide separation distance of 0.1 μm. This figure was generated using the [`plot2D`](../Python_User_Interface.md#data-visualization) routine and shows the source and flux monitor (red hatches), force monitors (blue lines), and PMLs (green hatches) surrounding the cell. From the force spectra shown above, at this separation distance the anti-symmetric mode is repulsive whereas the symmetric mode is attractive.
![](../images/parallel_wvgs_s0.1.png) diff --git a/doc/docs/Python_User_Interface.md b/doc/docs/Python_User_Interface.md index 016605ccd..1a97e847a 100644 --- a/doc/docs/Python_User_Interface.md +++ b/doc/docs/Python_User_Interface.md @@ -813,7 +813,7 @@ The index *n* (1,2,3,...) of the desired band ω*n*(**k**) to compute **`direction` [`mp.X`, `mp.Y`, or `mp.Z;` default `mp.AUTOMATIC`], `eig_match_freq` [`boolean;` default `True`], `eig_kpoint` [`Vector3`]** — -By default (if `eig_match_freq` is `True`), Meep tries to find a mode with the same frequency ω*n*(**k**) as the `src` property (above), by scanning **k** vectors in the given `direction` using MPB's `find_k` functionality. Alternatively, if `eig_kpoint` is supplied, it is used as an initial guess for **k**. By default, `direction` is the direction normal to the source region, assuming `size` is $d$–1 dimensional in a $d$-dimensional simulation (e.g. a plane in 3d). If `direction` is set to `mp.NO_DIRECTION`, then `eig_kpoint` 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 [launch modes in oblique ridge waveguides](Python_Tutorials/Eigenmode_Source.md#index-guided-modes-in-a-ridge-waveguide) (not perpendicular to the source plane). If `eig_match_freq` is `False`, then the specific **k** vector of the desired mode is specified with `eig_kpoint` (in Meep units of 2π/(unit length)). By default, the **k** components in the plane of the source region are zero. However, if the source 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 `eig_kpoint`. Note that once **k** is either found by MPB, or specified by `eig_kpoint`, the field profile used to create the current sources corresponds to the [Bloch mode](https://en.wikipedia.org/wiki/Bloch_wave), $\mathbf{u}_{n,\mathbf{k}}(\mathbf{r})$, multiplied by the appropriate exponential factor, $e^{i \mathbf{k} \cdot \mathbf{r}}$. +By default (if `eig_match_freq` is `True`), Meep tries to find a mode with the same frequency ω*n*(**k**) as the `src` property (above), by scanning **k** vectors in the given `direction` using MPB's `find_k` functionality. Alternatively, if `eig_kpoint` is supplied, it is used as an initial guess for **k**. By default, `direction` is the direction normal to the source region, assuming `size` is $d$–1 dimensional in a $d$-dimensional simulation (e.g. a plane in 3d). If `direction` is set to `mp.NO_DIRECTION`, then `eig_kpoint` 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 [launch modes in oblique ridge waveguides](Python_Tutorials/Eigenmode_Source.md#oblique-waveguides) (not perpendicular to the source plane). If `eig_match_freq` is `False`, then the **k** vector of the desired mode is specified with `eig_kpoint` (in Meep units of 2π/(unit length)). Also, the eigenmode frequency computed by MPB overwrites the `frequency` parameter of the `src` property for a `GaussianSource` and `ContinuousSource` but not `CustomSource` (the `width` or any other parameter of `src` is unchanged). By default, the **k** components in the plane of the source region are zero. However, if the source 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 `eig_kpoint`. Note that once **k** is either found by MPB, or specified by `eig_kpoint`, the field profile used to create the current sources corresponds to the [Bloch mode](https://en.wikipedia.org/wiki/Bloch_wave), $\mathbf{u}_{n,\mathbf{k}}(\mathbf{r})$, multiplied by the appropriate exponential factor, $e^{i \mathbf{k} \cdot \mathbf{r}}$. **`eig_parity` [`mp.NO_PARITY` (default), `mp.EVEN_Z`, `mp.ODD_Z`, `mp.EVEN_Y`, `mp.ODD_Y`]** — @@ -1254,6 +1254,8 @@ The parameters of this routine are the same as that of `get_eigenmode_coefficien + `kdom`: the dominant planewave of mode `band_num` + `amplitude(point, component)`: the (complex) value of the given E or H field `component` (`Ex`, `Hy`, etcetera) at a particular `point` (a `Vector3`) in space (interpreted with Bloch-periodic boundary conditions if you give a point outside the original `eig_vol`). +If `match_frequency=False` or `kpoint` is not zero in the given `direction`, the `frequency` input parameter is ignored. + **`get_eigenmode_freqs(flux)`** — Given a flux object, returns a list of the frequencies that it is computing the spectrum for. diff --git a/doc/docs/Scheme_Tutorials/Eigenmode_Source.md b/doc/docs/Scheme_Tutorials/Eigenmode_Source.md index c6d9a3cb1..f5505c409 100644 --- a/doc/docs/Scheme_Tutorials/Eigenmode_Source.md +++ b/doc/docs/Scheme_Tutorials/Eigenmode_Source.md @@ -89,7 +89,9 @@ For the multi-mode case, a constant-amplitude current source excites a superposi Note that in `eigenmode-source`, the `direction` property must be set to `NO-DIRECTION` for a non-zero `eig-kpoint` which specifies the waveguide axis. -Additionally, we can demonstrate the eigenmode source for a rotated waveguide. The results are shown in the two figures below for the single- and multi-mode case. There is one subtlety: for mode **A** in the multi-mode case, the `bnum` parameter is set to 3 rather than 2. This is because a non-zero rotation angle breaks the symmetry in the $y$-direction which therefore precludes the use of `EVEN-Y` in `eig-parity`. Without any parity specified for the $y$-direction, the second band corresponds to *odd* modes. This is why we must select the third band which contains even modes. +### Oblique Waveguides + +The eigenmode source can also be used to launch modes in an oblique/rotated waveguide. The results are shown in the two figures below for the single- and multi-mode case. There is one subtlety: for mode **A** in the multi-mode case, the `bnum` parameter is set to 3 rather than 2. This is because a non-zero rotation angle breaks the symmetry in the $y$-direction which therefore precludes the use of `EVEN-Y` in `eig-parity`. Without any parity specified for the $y$-direction, the second band corresponds to *odd* modes. This is why we must select the third band which contains even modes. Note that an oblique waveguide leads to a breakdown in the [PML](../Perfectly_Matched_Layer.md#breakdown-of-pml-in-inhomogeneous-media). A simple workaround for mitigating the PML reflection artifacts in this case is to double the `thickness` from 1 to 2. diff --git a/doc/docs/Scheme_Tutorials/Optical_Forces.md b/doc/docs/Scheme_Tutorials/Optical_Forces.md index 98977a857..2dd1529e7 100644 --- a/doc/docs/Scheme_Tutorials/Optical_Forces.md +++ b/doc/docs/Scheme_Tutorials/Optical_Forces.md @@ -10,7 +10,7 @@ This tutorial demonstrates Meep's ability to compute classical forces via the [M The gradient force $F$ on each waveguide arising from the evanescent coupling of the two waveguide modes can be computed analytically: -$$F=-\frac{1}{\omega}\frac{d\omega}{ds}\Bigg\vert_\vec{k}U,$$ +$$F=-\frac{1}{\omega}\frac{d\omega}{ds}\Bigg\vert_\vec{k}U$$ where $\omega$ is the mode frequency of the coupled waveguide system, $s$ is the separation distance between the parallel waveguides, $k$ is the conserved wave vector, and $U$ is the total energy of the electromagnetic fields. By convention, negative and positive values correspond to attractive and repulsive forces, respectively. For more details, see [Optics Letters, Vol. 30, pp. 3042-4, 2005](https://www.osapublishing.org/ol/abstract.cfm?uri=ol-30-22-3042). This expression has been shown to be mathematically equivalent to the Maxwell stress tensor in [Optics Express, Vol. 17, pp. 18116-35, 2009](http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-17-20-18116). We will verify this result in this tutorial. In this particular example, only the fundamental mode with odd mirror-symmetry in $y$ shows the bidirectional force. @@ -77,7 +77,7 @@ Since we do not know apriori what the mode frequency is for a given waveguide se (print "freq:, " s ", " f "\n") ``` -Once we have determined the mode frequency, we then replace the `source` with [`eigenmode-source`](../Scheme_User_Interface.md#eigenmode-source) to perform the main simulation: compute (1) the force on each waveguide due to the mode coupling and (2) the power in the mode. An `eigenmode-source` with `eig-match-freq?` set to `false` is then used to launch the guided mode. Alternatively, a constant-amplitude point/area source can be used to launch the mode but this is less efficient as demonstrated in [Tutorial/Eigenmode Source/Index-Guided Modes in a Ridge Waveguide](Eigenmode_Source.md#index-guided-modes-in-a-ridge-waveguide). +Once we have determined the mode frequency, we then replace the `source` with [`eigenmode-source`](../Scheme_User_Interface.md#eigenmode-source) to perform the main simulation: compute (1) the force on each waveguide due to the mode coupling and (2) the power in the mode. An `eigenmode-source` with `eig-match-freq?` set to `false` is then used to launch the guided mode using a pulse (the `frequency` but not the `fwidth` parameter of `gaussian-src` is ignored). Alternatively, a constant-amplitude point/area source can be used to launch the mode but this is less efficient as demonstrated in [Tutorial/Eigenmode Source/Index-Guided Modes in a Ridge Waveguide](Eigenmode_Source.md#index-guided-modes-in-a-ridge-waveguide). ```scm (reset-meep) @@ -113,7 +113,7 @@ In this example, the fields of the guided mode never decay away to zero. [Choosi The simulation is run over the range of separation distances $s$ from 0.05 to 1.00 μm in increments of 0.05 μm. The results are compared with those from MPB. This is shown in the top figure. The two methods show good agreement. -The following figure shows the $E_y$ mode profile at a waveguide separation distance of 0.1 μm. This figure shows the source/flux region (red hatches), force monitors (blue lines), and PMLs (green hatches) surrounding the cell. From the force spectra shown above, at this separation distance the anti-symmetric mode is repulsive whereas the symmetric mode is attractive. +The following figure shows the $E_y$ mode profiles at a waveguide separation distance of 0.1 μm. This figure shows the source and flux monitor (red hatches), force monitors (blue lines), and PMLs (green hatches) surrounding the cell. From the force spectra shown above, at this separation distance the anti-symmetric mode is repulsive whereas the symmetric mode is attractive.
![](../images/parallel_wvgs_s0.1.png) diff --git a/doc/docs/Scheme_User_Interface.md b/doc/docs/Scheme_User_Interface.md index 434a19690..aeeec4594 100644 --- a/doc/docs/Scheme_User_Interface.md +++ b/doc/docs/Scheme_User_Interface.md @@ -665,7 +665,7 @@ The index *n* (1,2,3,...) of the desired band ω*n*(**k**) to compute **`direction` [`X`, `Y`, or `Z;` default `AUTOMATIC`], `eig-match-freq?` [`boolean;` default `true`], `eig-kpoint` [`vector3`]** — -By default (if `eig-match-freq?` is `true`), Meep tries to find a mode with the same frequency ω*n*(**k**) as the `src` property (above), by scanning **k** vectors in the given `direction` using MPB's `find-k` functionality. Alternatively, if `eig-kpoint` is supplied, it is used as an initial guess for **k**. By default, `direction` is the direction normal to the source region, assuming `size` is $d$–1 dimensional in a $d$-dimensional simulation (e.g. a plane in 3d). If `direction` is set to `NO-DIRECTION`, then `eig_kpoint` 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 [launch modes in oblique ridge waveguides](Scheme_Tutorials/Eigenmode_Source.md#index-guided-modes-in-a-ridge-waveguide) (not perpendicular to the source plane). If `eig-match-freq?` is `false`, then the specific **k** vector of the desired mode is specified with `eig-kpoint` (in Meep units of 2π/(unit length)). By default, the **k** components in the plane of the source region are zero. However, if the source 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 `eig-kpoint`. Note that once **k** is either found by MPB, or specified by `eig-kpoint`, the field profile used to create the current sources corresponds to the [Bloch mode](https://en.wikipedia.org/wiki/Bloch_wave), $\mathbf{u}_{n,\mathbf{k}}(\mathbf{r})$, multiplied by the appropriate exponential factor, $e^{i \mathbf{k} \cdot \mathbf{r}}$. +By default (if `eig-match-freq?` is `true`), Meep tries to find a mode with the same frequency ω*n*(**k**) as the `src` property (above), by scanning **k** vectors in the given `direction` using MPB's `find-k` functionality. Alternatively, if `eig-kpoint` is supplied, it is used as an initial guess for **k**. By default, `direction` is the direction normal to the source region, assuming `size` is $d$–1 dimensional in a $d$-dimensional simulation (e.g. a plane in 3d). If `direction` is set to `NO-DIRECTION`, then `eig-kpoint` 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 [launch modes in oblique ridge waveguides](Scheme_Tutorials/Eigenmode_Source.md#oblique-waveguides) (not perpendicular to the source plane). If `eig-match-freq?` is `false`, then the **k** vector of the desired mode is specified with `eig-kpoint` (in Meep units of 2π/(unit length)). Also, the eigenmode frequency computed by MPB overwrites the `frequency` parameter of the `src` property for a `gaussian-src` and `continuous-src` but not `custom-src` (the `width` or any other parameter of `src` is unchanged). By default, the **k** components in the plane of the source region are zero. However, if the source 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 `eig-kpoint`. Note that once **k** is either found by MPB, or specified by `eig-kpoint`, the field profile used to create the current sources corresponds to the [Bloch mode](https://en.wikipedia.org/wiki/Bloch_wave), $\mathbf{u}_{n,\mathbf{k}}(\mathbf{r})$, multiplied by the appropriate exponential factor, $e^{i \mathbf{k} \cdot \mathbf{r}}$. **`eig-parity` [`NO-PARITY` (default), `EVEN-Z`, `ODD-Z`, `EVEN-Y`, `ODD-Y`]** — diff --git a/doc/docs/Subpixel_Smoothing.md b/doc/docs/Subpixel_Smoothing.md index ddc757021..d197df303 100644 --- a/doc/docs/Subpixel_Smoothing.md +++ b/doc/docs/Subpixel_Smoothing.md @@ -133,24 +133,24 @@ The adaptive numerical integration used for subpixel smoothing of material funct What Happens When Subpixel Smoothing is Disabled? ------------------------------------------------- -When subpixel smoothing is disabled by either (1) setting `eps_averaging=False` in the [`Simulation`](Python_User_Interface.md#the-simulation-class) constructor or (2) using a [material_function](Subpixel_Smoothing.md#enabling-averaging-for-material-function) (as is typical in the [adjoint solver](Python_Tutorials/AdjointSolver.md)), each electric field component ($E_x$, $E_y$, $E_z$) in a given voxel is individually assigned a scalar permittivity (for isotropic materials) based on whatever the value of the permittivity is at that position in the [Yee grid](Yee_Lattice.md). This results in [staircasing artifacts](Subpixel_Smoothing.md) due to the discontinuous material interfaces as well as the staggered nature of the Yee grid points. Any change in the resolution which shifts the location of the Yee grid points relative to the material interfaces will result in unpredictable changes to any computed quantities. The coordinates the Yee grid points can be obtained using a [field function](Field_Functions.md#coordinates-of-the-yee-grid) which can be useful for debugging. +When subpixel smoothing is disabled by either (1) setting `eps_averaging=False` in the [`Simulation`](Python_User_Interface.md#the-simulation-class) constructor or (2) using a [material function](Subpixel_Smoothing.md#enabling-averaging-for-material-function) (as is typical in the [adjoint solver](Python_Tutorials/AdjointSolver.md)), each electric field component ($E_x$, $E_y$, $E_z$) in a given voxel is individually assigned a scalar permittivity (for isotropic materials) based on whatever the value of the permittivity is at that position in the [Yee grid](Yee_Lattice.md). This results in [staircasing artifacts](Subpixel_Smoothing.md) due to the discontinuous material interfaces as well as the staggered nature of the Yee grid points. Any change in the resolution which shifts the location of the Yee grid points relative to the material interfaces will result in unpredictable changes to any computed quantities. The coordinates the Yee grid points can be obtained using a [field function](Field_Functions.md#coordinates-of-the-yee-grid) which can be useful for debugging. Subpixel Smoothing vs. Bilinear Interpolation --------------------------------------------- In certain cases, using subpixel smoothing may be impractical given the poor runtime performance of the adaptive numerical integration method as discussed previously. A partial workaround, to ensure that Meep responds continuously to changes in the simulation parameters (even if absolute accuracy is not improved) is to *interpolate* any discontinuous structure onto the Yee grid. Otherwise, tiny changes in Meep's [Yee grid](Yee_Lattice.md) due to e.g. small changes in the resolution could cause discontinuous jumps in ε. -As a demonstration of this effect, consider a ring resonator (same structure as [above](#continuously-varying-shapes-and-results)) in which the ring geometry can be represented using five different methods: +As a demonstration of this effect, consider a ring resonator (inner radius: 2 μm, width: 1 μm; same structure as [above](#continuously-varying-shapes-and-results)) in which the ring geometry can be represented using five different methods: -1. Two overlapping [`Cylinder`](Python_User_Interface.md#cylinder) objects (subpixel smoothing). -2. Two overlapping [`Prism`](Python_User_Interface.md#prism) objects (subpixel smoothing). Could also be imported as a [GDSII file](Python_User_Interface.md#gdsii-support). +1. two overlapping [`Cylinder`](Python_User_Interface.md#cylinder) objects (subpixel smoothing). +2. two overlapping [`Prism`](Python_User_Interface.md#prism) objects (subpixel smoothing). Could also be imported as a [GDSII file](Python_User_Interface.md#gdsii-support). 3. [material function](Python_User_Interface.md#material-function) (no smoothing). -4. Grid via [`epsilon_input_file`](Python_User_Interface.md#the-simulation-class) (no smoothing; bilinear interpolation onto Yee grid). -5. Grid via `epsilon_input_file` (no smoothing; no interpolation). +4. pixel grid via [`epsilon_input_file`](Python_User_Interface.md#the-simulation-class) (no smoothing; bilinear interpolation onto Yee grid). +5. pixel grid via `epsilon_input_file` (no smoothing; no interpolation). Of these five methods, (3) and (5) produce discontinuous structures. -The grid imported from `epsilon_input_file` via an HDF5 file in (4) and (5) is generated by the function `output_epsilon` when using the material function from (3) at a resolution of 80. +The pixel grid imported from the HDF5 `epsilon_input_file` in (4) and (5) is generated by the function `output_epsilon` when using the material function from (3) at a resolution of 80. The prism ring geometry in (2) is generated using: @@ -172,10 +172,10 @@ The following convergence plot shows the frequency for the resonant mode with $H ![](images/ring_freq_vs_resolution.png)
-There are three important items to note. (1) The grid and prism representations are each converging to a different frequency than the material function and cylinder. This is because in the limit of infinite resolution, they are *different* structures than the cylinders. (2) The material function is the same structure as the cylinder with no smoothing. In the limit of infinite resolution, the material function and cylinder converge to the same frequency. The only difference is the *rate* of convergence: the cylinder is second order (due to subpixel smoothing) whereas the material function is first order. See the convergence plot in the third figure from the top. (3) The non-interpolated grid shows irregular convergence compared with the interpolated grid. This is expected because the non-interpolated grid is discontinuous but the interpolated grid is not. Also, because these are different structures the two grids converge to different frequencies. To see this trend clearly requires reducing the "jumpiness" of the non-interpolated grid: the Meep resolution needs to be increased beyond 200 which is already ~3X the grid resolution. +There are three important items to note. (1) The pixel grid and prism representations are each converging to a different frequency than the material function and cylinder. This is because in the limit of infinite resolution, they are *different* structures than the cylinders. (2) The material function is the same structure as the cylinder with no smoothing. In the limit of infinite resolution, the material function and cylinder converge to the same frequency. The only difference is the *rate* of convergence: the cylinder is second order (due to subpixel smoothing) whereas the material function is first order. See the convergence plot above (third figure from the top). (3) The non-interpolated pixel grid shows irregular convergence compared with the interpolated grid. This is expected because the non-interpolated grid is discontinuous but the interpolated grid is not. Also, because these are different structures the two pixel grids converge to different frequencies. To see this trend clearly requires reducing the "jumpiness" of the non-interpolated grid: the Meep resolution needs to be increased beyond 200 which is already ~3X the grid resolution. -Since the interpolated grid has already been smoothed to a *continuous* $\varepsilon(x,y)$ function, subpixel smoothing (which is not supported for `epsilon_input_file`) is not really necessary once the Yee grid resolution exceeds the input image resolution. This can be seen in the above plot: for Meep Yee grid resolutions of 80 (equal to the grid resolution of the HDF5 file) and above, the changes in the results are much smaller than those at lower resolutions. Also, higher-order interpolation schemes are not necessary because the Yee discretization is already essentially equivalent to linear interpolation. +Since the interpolated pixel grid has already been smoothed to a continuous $\varepsilon(x,y)$ function, subpixel smoothing (which is not supported for `epsilon_input_file`) is not really necessary once the Yee grid resolution exceeds the input image resolution. This can be seen in the above plot: for Meep Yee grid resolutions of 80 (equal to the pixel grid resolution of the HDF5 file) and above, the changes in the results are much smaller than those at lower resolutions. Also, higher-order interpolation schemes are not necessary because the Yee discretization is already essentially equivalent to linear interpolation. -As a practical matter, increasing the Meep resolution beyond the resolution of a non-interpolated grid is not physically meaningful because this is trying to resolve the individual pixels of an imported image. In the case of a grid imported via `epsilon_input_file`, this is not an issue because the bilinear interpolation is performed automatically by default. However, no built-in interpolation is provided for a material function; it must be provided by the user (i.e., convolving the discontinuous material function with a smoothing kernel). As a corollary, when designing structures using a pixel grid (as in the [adjoint solver](Python_Tutorials/AdjointSolver.md)), the pixel density of the degrees of freedom should typically be at least as big as the Meep resolution if not greater. +As a practical matter, increasing the Meep resolution beyond the resolution of a non-interpolated pixel grid is not physically meaningful because this is trying to resolve the individual pixels of an imported image. In the case of a pixel grid imported via `epsilon_input_file`, this is not an issue because the bilinear interpolation is performed automatically by default. However, no built-in interpolation is provided for a material function; it must be provided by the user (i.e., convolving the discontinuous material function with a smoothing kernel). As a corollary, when designing structures using a pixel grid (as in the [adjoint solver](Python_Tutorials/AdjointSolver.md)), the pixel density of the degrees of freedom should typically be at least as big as the Meep resolution if not greater. In general, by making a discontinuous structure continuous, via subpixel smoothing or some other form of interpolation, the convergence becomes more regular (the results change more continuously to changes in resolution or other parameters), although it does not necessarily become more accurate compared to the desired infinite-resolution structure unless the full anisotropic smoothing is performed. If the initial structure is already continuous, no additional preprocessing is necessary. diff --git a/doc/docs/images/ring_freq_vs_resolution.png b/doc/docs/images/ring_freq_vs_resolution.png index de6cd3d32..02bc23d09 100644 Binary files a/doc/docs/images/ring_freq_vs_resolution.png and b/doc/docs/images/ring_freq_vs_resolution.png differ diff --git a/python/examples/parallel-wvgs-force.ipynb b/python/examples/parallel-wvgs-force.ipynb index b3bdbff88..32e9eb0d8 100644 --- a/python/examples/parallel-wvgs-force.ipynb +++ b/python/examples/parallel-wvgs-force.ipynb @@ -12,7 +12,7 @@ "\n", "The gradient force $F$ on each waveguide arising from the evanescent coupling of the two waveguide modes can be computed analytically:\n", "\n", - "$$F=-\\frac{1}{\\omega}\\frac{d\\omega}{ds}\\Bigg\\vert_\\vec{k}U,$$\n", + "$$F=-\\frac{1}{\\omega}\\frac{d\\omega}{ds}\\Bigg\\vert_\\vec{k}U$$\n", "\n", "where $\\omega$ is the mode frequency of the coupled waveguide system, $s$ is the separation distance between the parallel waveguides, $k$ is the conserved wave vector, and $U$ is the total energy of the electromagnetic fields. By convention, negative and positive values correspond to attractive and repulsive forces, respectively. For more details, see [Optics Letters, Vol. 30, pp. 3042-4, 2005](https://www.osapublishing.org/ol/abstract.cfm?uri=ol-30-22-3042). This expression has been shown to be mathematically equivalent to the Maxwell stress tensor in [Optics Express, Vol. 17, pp. 18116-35, 2009](http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-17-20-18116). We will verify this result in this tutorial. In this particular example, only the fundamental mode with odd mirror-symmetry in $y$ shows the bidirectional force.\n", "\n", @@ -20,7 +20,7 @@ "\n", "The gradient force $F$ can be computed using two different methods: (1) using MPB, compute the frequency $\\omega$ and group velocity $v_g$ for a given mode over a range of separation distances $s$ and then use a centered [finite-difference](https://en.wikipedia.org/wiki/Finite_difference) scheme to evaluate $F$ using the formula from above, and (2) using Meep, directly compute both the gradient force $F$ and the transmitted power $P$ over the same range of separation distances $s$. This tutorial verifies that (1) and (2) produce equivalent results.\n", "\n", - "The main component of the Meep script is the function `parallel_waveguide` which computes the force $F$ and transmitted power $P$ given the waveguide separation distance `s` and relative phase of the waveguide modes `xodd=True/False`. The eigenmode frequency $\\omega$ is computed first using `get_eigenmode` in order to specify the frequency of the `add_force` and `add_flux` monitors. An [`EigenModeSource`](https://meep.readthedocs.io/en/latest/Python_User_Interface/#eigenmodesource) with `eig_match_freq=False` is then used to launch the guided mode. Alternatively, a constant-amplitude point/area source can be used to launch the mode but this is less efficient as demonstrated in [Tutorial/Eigenmode Source/Index-Guided Modes in a Ridge Waveguide](https://meep.readthedocs.io/en/latest/Python_Tutorials/Eigenmode_Source/#index-guided-modes-in-a-ridge-waveguide). The waveguide has width/height of $a=1$ μm and a fixed propagation wavevector of $\\pi/a$." + "The main component of the Meep script is the function `parallel_waveguide` which computes the force $F$ and transmitted power $P$ given the waveguide separation distance `s` and relative phase of the waveguide modes `xodd=True/False`. The eigenmode frequency $\\omega$ is computed first using `get_eigenmode` in order to specify the frequency of the `add_force` and `add_flux` monitors. (Note that when `match_frequency=False`, `get_eigenmode` ignores the input `frequency` parameter.) An [`EigenModeSource`](https://meep.readthedocs.io/en/latest/Python_User_Interface/#eigenmodesource) with `eig_match_freq=False` is then used to launch the guided mode using a pulse (the `frequency` but not the `fwidth` parameter of `GaussianSource` is ignored). Alternatively, a constant-amplitude point/area source can be used to launch the mode but this is less efficient as demonstrated in [Tutorial/Eigenmode Source/Index-Guided Modes in a Ridge Waveguide](https://meep.readthedocs.io/en/latest/Python_Tutorials/Eigenmode_Source/#index-guided-modes-in-a-ridge-waveguide). The waveguide has width/height of $a=1$ μm and a fixed propagation wavevector of $\\pi/a$." ] }, { @@ -2363,11 +2363,9 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "The following figure shows the $E_y$ mode profile at a waveguide separation distance of 0.1 μm. This figure was generated using the [`plot2D`](https://meep.readthedocs.io/en/latest/Python_User_Interface/#data-visualization) routine and shows the source/flux region (red hatches), force monitors (blue lines), and PMLs (green hatches) surrounding the cell. From the force spectra shown above, at this separation distance the anti-symmetric mode is repulsive whereas the symmetric mode is attractive.\n", + "The following figure shows the $E_y$ mode profiles at a waveguide separation distance of 0.1 μm. This figure was generated using the [`plot2D`](https://meep.readthedocs.io/en/latest/Python_User_Interface/#data-visualization) routine and shows the source and flux monitor (red hatches), force monitors (blue lines), and PMLs (green hatches) surrounding the cell. From the force spectra shown above, at this separation distance the anti-symmetric mode is repulsive whereas the symmetric mode is attractive.\n", "\n", - "
\n", "![](https://meep.readthedocs.io/en/latest/images/parallel_wvgs_s0.1.png)\n", - "
\n", "\n", "The MPB simulation is in [examples/parallel-wvgs-mpb.py](https://github.com/NanoComp/meep/blob/master/python/examples/parallel-wvgs-mpb.py). There are important differences related to the coordinate dimensions between the MPB and Meep scripts. In the MPB script, the 2d cell is defined using the $yz$ plane, the waveguide propagation axis is $x$,and the waveguide separation axis is $y$. As a consequence, the `num_bands` parameter is always `1` since the $y$ parity of the mode can be defined explicitly (i.e., `run_yodd_zodd` vs. `run_yeven_zodd`). This is different from the Meep script since Meep requires that a 2d cell be defined in the $xy$ plane. MPB has no such requirement." ] diff --git a/src/mpb.cpp b/src/mpb.cpp index cd507e255..66ada5ab8 100644 --- a/src/mpb.cpp +++ b/src/mpb.cpp @@ -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(frequency); + if (!match_frequency) src_mpb->set_frequency(global_eigenmode_data->frequency); /*--------------------------------------------------------------*/ // step 2: add sources whose radiated field reproduces the */