Meep supports the simulation of Maxwell's equations in cylindrical coordinates for structures that have continuous rotational symmetry around the z axis. This reduces problems in 3d to 2d, and 2d to 1d, if there is sufficient symmetry.
[TOC]
In Tutorial/Basics/Modes of a Ring Resonator, the modes of a ring resonator were computed by performing a 2d simulation. This example involves simulating the same structure while exploiting the fact that the system has continuous rotational symmetry, by performing the simulation in cylindrical coordinates. The simulation script is in examples/ring-cyl.ctl.
The parameters of the problem are defined with exactly the same values as in the 2d simulation:
(define-param n 3.4) ; index of waveguide
(define-param w 1) ; width of waveguide
(define-param r 1) ; inner radius of ring
(define-param pad 4) ; padding between waveguide and edge of PML
(define-param dpml 2) ; thickness of PML
The dimensions and size of the computational cell are defined:
(define sr (+ r w pad dpml)) ; radial size (cell is from 0 to sr)
(set! dimensions CYLINDRICAL)
(set! geometry-lattice (make lattice (size sr no-size no-size)))
The key thing is to set the dimensions
parameter to CYLINDRICAL
. This means that all vectors represent (sr = r + w + pad + dpml
, and runs from 0
to sr
(by default) rather than from -sr/2
to sr/2
as it would for any other dimension. Note that the no-size
because it is in 2d. The no-size
, corresponding to the continuous rotational symmetry. A finite
In particular, in systems with continuous rotational symmetry, by an analogue of Bloch's theorem, the angular dependence of the fields can always be chosen in the form m
which is set to a command-line argument that is 3 by default.
(set-param! m 3)
This is essentially a 1d calculation where Meep must discretize the
The geometry is now specified by a single block
object — remember that this is a block in cylindrical coordinates, so that it really specifies an annular ring:
(set! geometry (list
(make block (center (+ r (/ w 2))) (size w infinity infinity)
(material (make dielectric (index n))))))
(set! pml-layers (list (make pml (thickness dpml))))
(set-param! resolution 10)
PMLs are on "all" sides. The
The remaining inputs are almost exactly the same as in the previous 2d simulation. A single Gaussian point source is added in the
(define-param fcen 0.15) ; pulse center frequency
(define-param df 0.1) ; pulse width (in frequency)
(set! sources (list
(make source
(src (make gaussian-src (frequency fcen) (fwidth df)))
(component Ez) (center (+ r 0.1) 0 0))))
Note that this isn't really a point source, however, because of the cylindrical symmetry — it is really a ring source with
(run-sources+ 200 (after-sources (harminv Ez (vector3 (+ r 0.1)) fcen df)))
At the very end, one period of the fields is output to create an animation. A single field output would be a 1d dataset along the to-appended
is used to append these datasets to a single HDF5 file to obtain an in-volume
is used to specify a larger output volume than just the computational cell: in particular, the output is from -sr
to sr
in the
(run-until (/ 1 fcen)
(in-volume (volume (center 0) (size (* 2 sr)))
(at-beginning output-epsilon)
(to-appended "ez"
(at-every (/ 1 fcen 20) output-efield-z))))
The simulation is ready to be run. Recall that, in the 2d calculation, three modes were obtained in this frequency range: (1)
In particular, the three calculations are:
unix% meep m=3 ring-cyl.ctl | grep harminv
unix% meep m=4 ring-cyl.ctl | grep harminv
unix% meep m=5 ring-cyl.ctl | grep harminv
giving the combined output:
harminv0:, frequency, imag. freq., Q, |amp|, amplitude, error
harminv0:, 0.11835455441250631, -6.907792691647415e-4, 85.66741917111612, 0.02570190626349302, -0.02402703883357199-0.00912630212448642i, 5.286949731053267e-10+0.0i
harminv0:, 0.1475578747705309, -1.938438860632441e-4, 380.61008208014414, 0.19361245519715206, 0.1447225471614173+0.12861246887677943i, 5.889273063545974e-11+0.0i
harminv0:, 0.1759448592380757, -4.900590034953583e-5, 1795.1395442502285, 0.0452479314013276, -0.014395016792255884-0.042897072017212545i, 1.6343462235932872e-10+0.0i
This is indeed very close to the 2d simulations: the frequencies are within 1% of the previous values. The
Which is more accurate, the 2d or the cylindrical simulation? This question can be answered by increasing the resolution in both cases and seeing what they converge towards. In particular, let's focus on the
The problem is twofold. First, there is some signal-processing error in determining harminv
output which is only 4e-7 for the 2d simulation vs. 6e-11 for the cylindrical case. This error can be reduced by running with a narrower bandwidth source, which excites just one mode and gives a cleaner signal, or by analyzing over a longer time than 200. Doing the former, we find that the 2d value of pad
) and/or increasing the PML thickness (increasing dpml
) so that it turns on more adiabatically. In the cylindrical simulation for resolution = 20
, if the PML thickness is increased to dpml = 16
, the dpml = 32
the
This illustrates the general principle that you need to check several parameters to ensure that results are converged in time-domain simulations: the resolution, the run time, the PML thickness, etcetera.
Finally, the field images are obtained. Since one mode per m
is being excited here anyway, according to harminv
, there is no real need for a narrow-band source. This will be used anyway just to remind you of the general procedure, however, e.g. for the
unix% meep m=3 fcen=0.118 df=0.01 ring-cyl.ctl
unix% h5topng -S 2 -Zc dkbluered -C ring-cyl-eps-001200.00.h5 ring-cyl-ez.h5
Note that, because of the to-appended
command, the ring-cyl-ez.h5
file is a 160$\times$18 dataset corresponding to an
Because only the
For a given mode of the ring resonator, it is often useful to know how sensitive the resonant frequency
Pertubation theory for Maxwell equations involving high index-contrast dielectric interfaces is reviewed in Chapter 2 of Photonics Crystals: Molding the Flow of Light, 2nd Edition (2008). The formula (equation 30 on p.19) for the frequency shift
$$ \Delta\omega = -\frac{\omega}{2} \frac{ \iint d^2 \vec{r} \big[ (\varepsilon_1 - \varepsilon_2) |\vec{E}{\parallel}(\vec{r})|^2 - \big(\frac{1}{\varepsilon_1} - \frac{1}{\varepsilon_2}\big)|\varepsilon\vec{E}{\perp}|^2\big] \Delta h}{\int d^3\vec{r} \varepsilon(\vec{r})|\vec{E}(\vec{r})|^2} + O(\Delta h^2)$$
In this expression, $\vec{E}{\parallel}(\vec{r})$ is the component of $\vec{E}$ that is parallel to the surface, and $\varepsilon\vec{E}{\perp}$ is the component of
There are three parts to the calculation: (1) find the resonant frequency of the unperturbed geometry using a broadband Gaussian source, (2) find the resonant mode profile of the unperturbed geometry using a narrowband source and from these fields compute the gradient via the perturbation-theory formula, and (3) find the resonant frequency of the perturbed geometry and from this compute the gradient using the finite-difference approximation. The perturbation is applied only to the inner and outer ring radii. The ring width is constant. There are two types of modes which are computed in separate simulations using different source polarizations: parallel (
The simulation script is in examples/perturbation-theory.ctl.
(set-param! resolution 100) ; pixels/um
(define-param perpendicular? true)
(define src-cmpt (if perpendicular? Hz Ez))
(define fcen (if perpendicular? 0.21 0.17)) ; pulse center frequency
(define-param n 3.4) ; index of waveguide
(define-param w 1) ; width of waveguide
(define-param r 1) ; inner radius of ring
(define-param pad 4) ; padding between waveguide and edge of PML
(define-param dpml 2) ; thickness of PML
(set-param! m 5) ; angular dependence
(set! pml-layers (list (make pml (thickness dpml))))
(define sr (+ r w pad dpml)) ; radial size (cell is from 0 to sr)
(set! dimensions CYLINDRICAL)
(set! geometry-lattice (make lattice (size sr no-size no-size)))
(set! geometry (list (make block
(center (+ r (/ w 2)))
(size w infinity infinity)
(material (make dielectric (index n))))))
;; find resonant frequency of unperturbed geometry using broadband source
(define df (* 0.2 fcen)) ; pulse width (in frequency)
(set! sources (list
(make source
(src (make gaussian-src (frequency fcen) (fwidth df)))
(component src-cmpt)
(center (+ r 0.1)))))
(run-sources+ 100 (after-sources (harminv src-cmpt (vector3 (+ r 0.1)) fcen df)))
(define frq-unperturbed (harminv-freq-re (car harminv-results)))
(reset-meep)
;; unperturbed geometry with narrowband source centered at resonant frequency
(set! pml-layers (list (make pml (thickness dpml))))
(set! geometry-lattice (make lattice (size sr no-size no-size)))
(set! geometry (list (make block
(center (+ r (/ w 2)))
(size w infinity infinity)
(material (make dielectric (index n))))))
(set! fcen frq-unperturbed)
(set! df (* 0.05 fcen))
(set! sources (list
(make source
(src (make gaussian-src (frequency fcen) (fwidth df)))
(component src-cmpt)
(center (+ r 0.1)))))
(run-sources+ 100)
(define deps (- 1 (* n n)))
(define deps-inv (- 1 (/ (* n n))))
(define numerator-integral 0)
(if perpendicular?
(let ((para-integral (* deps 2 pi (- (* r (sqr (magnitude (get-field-point Ep (vector3 r))))) (* (+ r w) (sqr (magnitude (get-field-point Ep (vector3 (+ r w)))))))))
(perp-integral (* deps-inv 2 pi (- (* (+ r w) (sqr (magnitude (get-field-point Dr (vector3 (+ r w)))))) (* r (sqr (magnitude (get-field-point Dr (vector3 r)))))))))
(set! numerator-integral (+ para-integral perp-integral)))
(set! numerator-integral (* deps 2 pi (- (* r (sqr (magnitude (get-field-point Ez (vector3 r))))) (* (+ r w) (sqr (magnitude (get-field-point Ez (vector3 (+ r w))))))))))
(define denominator-integral (electric-energy-in-box (volume (center (* 0.5 sr)) (size sr))))
(define perturb-theory-dw-dR (/ (* -1 frq-unperturbed numerator-integral) (* 4 denominator-integral)))
(reset-meep)
;; perturbed geometry with narrowband source
(define-param dr 0.01)
(set! pml-layers (list (make pml (thickness dpml))))
(set! geometry-lattice (make lattice (size sr no-size no-size)))
(set! geometry (list (make block
(center (+ r dr (/ w 2)))
(size w infinity infinity)
(material (make dielectric (index n))))))
(set! sources (list
(make source
(src (make gaussian-src (frequency fcen) (fwidth df)))
(component src-cmpt)
(center (+ r dr 0.1)))))
(run-sources+ 100 (after-sources (harminv src-cmpt (vector3 (+ r 0.1)) fcen df)))
(define frq-perturbed (harminv-freq-re (car harminv-results)))
(define finite-diff-dw-dR (/ (- frq-perturbed frq-unperturbed) dr))
(print "dwdR:, " perturb-theory-dw-dR " (pert. theory), " finite-diff-dw-dR " (finite diff.)\n")
There are three things to note. First, there is a built-in function electric-energy-in-box
which calculates the integral of get-field-point
are automatically interpolated; i.e., no additional post-processing is necessary. Third, when comparing the perturbation-theory result to the finite-difference approximation, there are two convergence parameters: the resolution and
For an resolution = 100
the results are:
dwdR:, -0.08544696397218979 (pert. theory), -0.08521249090736038 (finite diff.)
Doubling the resolution to 200, the results are:
dwdR:, -0.08544607322081005 (pert. theory), -0.08521153501551137 (finite diff.)
Both results have converged to at least five digits. The relative error at resolution 200 is 0.3%. The mode has a
For an resolution = 100
the results are:
dwdR:, -0.0805038571770864 (pert. theory), -0.07980873307536773 (finite diff.)
Doubling the resolution to 200, the results are:
dwdR:, -0.08020283464036788 (pert. theory), -0.07980880151594316 (finite diff.)
Both results have converged to at least three digits. The relative error at resolution 200 is 0.5%. The error is larger in this case due to the presence of the discontinuous fields at the dielectric interface. The mode has a
Finally, as reference, the same calculation can be set up in Cartesian coordinates as a 2d simulation. The simulation script is in examples/perturbation-theory-2d.ctl. There is one major difference in the 2d calculation: the mode produced by a point source in 2d is actually the
As an alternative to the "ring" sources of the previous examples, it is also possible to launch planewaves in cylindrical coordinates. This is demonstrated in this example which involves computing the scattering cross section of a finite-height dielectric cylinder. The results for the 2d simulation involving the cylindrical (
The calculation of the scattering cross section is described in Tutorial/Basics/Mie Scattering of a Lossless Dielectric Sphere which is modified for this example. A linearly-polarized (
(Note: a
In practice, this involves performing two separate simulations for
The simulation script is in examples/cylinder-cross-section.ctl.
(define-param r 0.7) ;; radius of cylinder
(define-param h 2.3) ;; height of cylinder
(define wvl-min (/ (* 2 pi r) 10))
(define wvl-max (/ (* 2 pi r) 2))
(define frq-min (/ wvl-max))
(define frq-max (/ wvl-min))
(define frq-cen (* 0.5 (+ frq-min frq-max)))
(define dfrq (- frq-max frq-min))
(define nfrq 100)
;; at least 8 pixels per smallest wavelength, i.e. (floor (/ 8 wvl-min))
(set-param! resolution 25)
(define dpml (* 0.5 wvl-max))
(define dair (* 1.0 wvl-max))
(define boundary-layers (list (make pml (thickness dpml))))
(set! pml-layers boundary-layers)
(define sr (+ r dair dpml))
(define sz (+ dpml dair h dair dpml))
(define cell (make lattice (size sr 0 sz)))
(set! geometry-lattice cell)
(set! dimensions CYLINDRICAL)
(set-param! m -1)
;; (is-integrated? true) necessary for any planewave source extending into PML
(define circ-src (list (make source
(src (make gaussian-src (frequency frq-cen) (fwidth dfrq) (is-integrated? true)))
(center (* 0.5 sr) 0 (+ (* -0.5 sz) dpml))
(size sr 0 0)
(component Er))
(make source
(src (make gaussian-src (frequency frq-cen) (fwidth dfrq) (is-integrated? true)))
(center (* 0.5 sr) 0 (+ (* -0.5 sz) dpml))
(size sr 0 0)
(component Ep)
(amplitude 0-1i))))
(set! sources circ-src)
(define box-z1 (add-flux frq-cen dfrq nfrq
(make flux-region (center (* 0.5 r) 0 (* -0.5 h)) (size r 0 0))))
(define box-z2 (add-flux frq-cen dfrq nfrq
(make flux-region (center (* 0.5 r) 0 (* +0.5 h)) (size r 0 0))))
(define box-r (add-flux frq-cen dfrq nfrq
(make flux-region (center r 0 0) (size 0 0 h))))
(run-sources+ 10)
(display-fluxes box-z1)
(save-flux "box-z1-flux" box-z1)
(save-flux "box-z2-flux" box-z2)
(save-flux "box-r-flux" box-r)
(reset-meep)
(define ncyl 2.0)
(set! geometry (list
(make block
(material (make medium (index ncyl)))
(size r 0 h)
(center (* 0.5 r) 0 0))))
(set! geometry-lattice cell)
(set! pml-layers boundary-layers)
(set! sources circ-src)
(define box-z1 (add-flux frq-cen dfrq nfrq
(make flux-region (center (* 0.5 r) 0 (* -0.5 h)) (size r 0 0))))
(define box-z2 (add-flux frq-cen dfrq nfrq
(make flux-region (center (* 0.5 r) 0 (* +0.5 h)) (size r 0 0))))
(define box-r (add-flux frq-cen dfrq nfrq
(make flux-region (center r 0 0) (size 0 0 h))))
(load-minus-flux "box-z1-flux" box-z1)
(load-minus-flux "box-z2-flux" box-z2)
(load-minus-flux "box-r-flux" box-r)
(run-sources+ 100)
(display-fluxes box-z1 box-z2 box-r)
Note that the "closed" DFT flux box is comprised of just three flux objects: two along display-fluxes
which computes the integral of the Poynting vector does so over the annular volume in cylindrical coordinates. There is no need for additional post-processing of the flux values.
The following Bash shell script runs the simulation. The script pipes the output to a file and extracts the input and and scattering flux data to separate files.
#!/bin/bash
meep cylinder-cross-section.ctl |tee cyl.out
grep flux1: cyl.out |cut -d, -f2- > input.dat
grep flux2: cyl.out |cut -d, -f2- > scatt.dat
The scattering cross section is computed from the simulation data and plotted using the following Matlab/Octave script.
input = dlmread('input.dat',',');
scatt = dlmread('scatt.dat',',');
freqs = input(:,1);
scatt_flux = scatt(:,2) - scatt(:,3) - scatt(:,4);
r = 0.7;
intensity = input(:,2)/(pi*r^2);
scatt_cross_section = -scatt_flux./intensity;
loglog(2*pi*r*freqs,scatt_cross_section,'bo-');
xlabel('(cylinder circumference)/wavelength, 2πr/\lambda');
ylabel('scattering cross section, \sigma');
title("Scattering Cross Section of a Lossless Dielectric Cylinder");
axis tight;
set(gca, "xminorgrid", "on");
set(gca, "yminorgrid", "on");
As shown below, the results for the scattering cross section computed using cylindrical coordinates agree well with the 3d Cartesian simulation. However, there is a large discrepancy in performance: for a single Intel Xeon 4.2GHz processor, the runtime of the cylindrical simulation is nearly 90 times shorter than the 3d simulation.
![](../images/cylinder_cross_section.png)It is also possible to compute a near-to-far field transformation in cylindrical coordinates. This is demonstrated in this example for a binary-phase zone plate which is a rotationally-symmetric diffractive lens used to focus a normally-incident planewave to a single spot.
Using scalar theory, the radius of the
where
The simulation script is in examples/zone-plate.ctl.
(set-param! resolution 25) ;; pixels/um
(define-param dpml 1.0) ;; PML thickness
(define-param dsub 2.0) ;; substrate thickness
(define-param dpad 2.0) ;; padding between zone plate and PML
(define-param zh 0.5) ;; zone-plate height
(define-param zN 25) ;; number of zones (odd zones: pi phase shift, even zones: none)
(define-param focal-length 200) ;; focal length of zone plate
(define-param spot-length 100) ;; far-field line length
(define-param ff-res 10) ;; far-field resolution (points/um)
(set! pml-layers (list (make pml (thickness dpml))))
(define-param wvl-cen 0.5)
(define frq-cen (/ wvl-cen))
(define dfrq (* 0.2 frq-cen))
(define r (map (lambda (n)
(sqrt (* n wvl-cen (+ focal-length (/ (* n wvl-cen) 4)))))
(arith-sequence 1 1 zN)))
(define sr (+ (list-ref r (- zN 1)) dpad dpml))
(define sz (+ dpml dsub zh dpad dpml))
(define cell (make lattice (size sr 0 sz)))
(set! geometry-lattice cell)
(set! dimensions CYLINDRICAL)
(set-param! m -1)
(set! sources (list (make source
(src (make gaussian-src (frequency frq-cen) (fwidth dfrq) (is-integrated? true)))
(center (* 0.5 sr) 0 (+ (* -0.5 sz) dpml))
(size sr 0 0)
(component Er))
(make source
(src (make gaussian-src (frequency frq-cen) (fwidth dfrq) (is-integrated? true)))
(center (* 0.5 sr) 0 (+ (* -0.5 sz) dpml))
(size sr 0 0)
(component Ep)
(amplitude 0-1i))))
(define glass (make medium (index 1.5)))
(set! geometry (list (make block
(material glass)
(size sr 0 (+ dpml dsub))
(center (* 0.5 sr) 0 (+ (* -0.5 sz) (* 0.5 (+ dpml dsub)))))))
(set! geometry (append geometry
(map (lambda (n)
(make block
(material (if (= (modulo n 2) 0) glass vacuum))
(size (list-ref r n) 0 zh)
(center (* 0.5 (list-ref r n)) 0 (+ (* -0.5 sz) dpml dsub (* 0.5 zh)))))
(reverse (arith-sequence 0 1 zN)))))
;; near-field monitor
(define n2f-obj (add-near2far frq-cen 0 1
(make near2far-region (center (* 0.5 (- sr dpml)) 0 (- (* 0.5 sz) dpml)) (size (- sr dpml) 0 0))
(make near2far-region (center (- sr dpml) 0 (- (* 0.5 sz) (* 0.5 (+ dsub zh dpad)))) (size 0 0 (+ dsub zh dpad)))))
(run-sources+ 100)
(output-farfields n2f-obj "r" (volume (center (* 0.5 (- sr dpml)) 0 (+ (* -0.5 sz) dpml dsub zh focal-length)) (size (- sr dpml) 0 0)) ff-res)
(output-farfields n2f-obj "z" (volume (center 0 0 (+ (* -0.5 sz) dpml dsub zh focal-length)) (size 0 0 spot-length)) ff-res)
Note that the volume specified in output-farfields
via center
and size
is in cylindrical coordinates. These points must therefore lie in the output-farfields
can be thought of as either cylindrical (output-farfields
tends to gradually slow down as the far-field point gets closer to the near-field monitor. This performance degradation is unavoidable and is due to the larger number of
The energy density of the far fields is computed from the simulation data and plotted using the following Matlab/Octave script.
focal_length = 200;
spot_length = 100;
sr = 53.3891;
dpml = 1.0;
subplot(1,2,1);
load "zone-plate-r.h5";
E2_r = abs(ex_r+1j*ex_i).^2 + abs(ey_r+1j*ey_i).^2 + abs(ez_r+1j*ez_i).^2;
r = linspace(0,sr-dpml,length(E2_r));
semilogy(r,E2_r,'bo-');
xlim([-2 20]);
xlabel('r coordinate (um)');
ylabel('energy density of far fields, |E|^2');
set (gca, "yminorgrid", "on");
subplot(1,2,2);
load "zone-plate-z.h5";
E2_z = abs(ex_r+1j*ex_i).^2 + abs(ey_r+1j*ey_i).^2 + abs(ez_r+1j*ez_i).^2;
z = linspace(focal_length-0.5*spot_length,focal_length+0.5*spot_length,length(E2_z));
semilogy(z,E2_z,'bo-');
xlabel('z coordinate (um)');
ylabel('energy density of far fields, |E|^2');
set (gca, "yminorgrid", "on");
Shown below is the far-field energy-density profile around the focal length for both the r and z coordinate directions for three lens designs with