perturbation-theory tutorial example

[stevengj]

It might be nice to add a tutorial example to https://meep.readthedocs.io/en/latest/Subpixel_Smoothing/ showing how to get the derivative dω/dR via perturbation theory rather than by center-difference approximations.

Basically, you just do the perturbation-theory surface integral with get_field at N=2πR/resolution equally spaced points around the ring's inner and outer surfaces — the average (multiplied by 2πR) is a good approximation for the surface integral.

[smartalecH]

In the spirit of perturbation theory, It might also be beneficial to add a tutorial (in an entirely new section) that calculates bent waveguide modes (with cylindrical symmetry) by perturbing the corresponding straight waveguide mode, as described in your paper.

Or maybe a tutorial that calculates loss from the lossless mode via perturbation theory [e.g. Snyder, Love].

This seems like a useful (and relatively straightforward) extension of the MPB tutorials.

[danielwboyce]

I've been working on some code for this tutorial, but have a couple of questions:

• While I'm pretty sure the surface integral from the perturbation theory paper is Eq. (12), how do we get from Eq. (12) to Eq. (3), which calculates dω/dR? In other words how do we get from the fields to dω/dR?
• Also what is meant by N in the last bullet point? With the way it's defined N<1 in a lot of cases, so do you just have one point on each of the ring's surfaces to check for the field?

from __future__ import division

import meep as mp
import numpy as np
from math import ceil
from statistics import mean

def main():
    n = 3.4                 # index of waveguide
    r = 1
    a = r                   # inner radius of ring
    w = 1                   # width of waveguide
    b = a+w                 # outer radius of ring
    pad = 4                 # padding between waveguide and edge of PML
    dpml = 2                # thickness of PML
    pml_layers = [mp.PML(thickness=dpml)]

    sxy = 2*(b+pad+dpml)  # cell size
    cell_size = mp.Vector3(sxy,sxy)
    resolution = 10

    c1 = mp.Cylinder(radius=b, material=mp.Medium(index=n))
    c2 = mp.Cylinder(radius=a)
    geometry = [c1,c2]

    # Exciting the first resonance mode (calculated with Harminiv)

    fcen = 0.118101372125519    # pulse center frequency
    df = 0.1                    # pulse width (in frequency)

    sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df), mp.Ez, mp.Vector3(r+0.1))]

    symmetries = [mp.Mirror(mp.Y)]

    sim = mp.Simulation(cell_size=cell_size,
                        geometry=geometry,
                        sources=sources,
                        resolution=resolution,
                        symmetries=symmetries,
                        boundary_layers=pml_layers)

    sim.run(mp.at_beginning(mp.output_epsilon),
            mp.after_sources(mp.Harminv(mp.Ez, mp.Vector3(r+0.1), fcen, df)),
            until_after_sources=300)

    # npts_inner = ceil(2 * np.pi * a / resolution)
    # npts_outer = ceil(2 * np.pi * b / resolution)
    npts_inner = 10
    npts_outer = 10
    angles_inner = 2 * np.pi / npts_inner * np.arange(npts_inner)
    angles_outer = 2 * np.pi / npts_outer * np.arange(npts_outer)

    inner_ring_fields = []
    outer_ring_fields = []
    for angle in angles_inner:
        point = mp.Vector3(a * np.cos(angle), a * np.sin(angle), 0)
        # what time step are these get_field_point calls coming from?
        e_x_field = sim.get_field_point(mp.Ex, point)
        e_y_field = sim.get_field_point(mp.Ey, point)
        e_z_field = sim.get_field_point(mp.Ez, point)
        e_total_field = np.real(np.sqrt(e_x_field*np.conj(e_x_field) + e_y_field*np.conj(e_y_field) + e_z_field*np.conj(e_z_field)))
        #e_total_field = np.real(np.sqrt(e_z_field * np.conj(e_z_field)))
        inner_ring_fields.append(e_total_field)

    for angle in angles_outer:
        point = mp.Vector3(a * np.cos(angle), a * np.sin(angle), 0)
        e_x_field = sim.get_field_point(mp.Ex, point)
        e_y_field = sim.get_field_point(mp.Ey, point)
        e_z_field = sim.get_field_point(mp.Ez, point)
        e_total_field = np.real(np.sqrt(e_x_field*np.conj(e_x_field) + e_y_field*np.conj(e_y_field) + e_z_field*np.conj(e_z_field)))
        #e_total_field = np.real(np.sqrt(e_z_field * np.conj(e_z_field)))
        outer_ring_fields.append(e_total_field)

    surface_integral = 2 * np.pi * b * (mean(inner_ring_fields) + mean(outer_ring_fields)) / 2
    print(f'\nThe value of the surface_integral is {surface_integral}')

if __name__ == '__main__':
    main()

[stevengj]

Note that the perturbation-theory formula is presented more simply in our book, in equation (30) on page 19.

[stevengj]

As explained in this tutorial, once you find the frequencies with Harminv, you run it again with a narrow-bandwidth source to excite just that mode. Then at the end of the second simulation you have the fields.

See the tutorial at "To see the mode, we will simply run the simulation again with a narrow-band source, …"

[stevengj]

In the manual you will see a function electric_energy_in_box. This calculates the integral of E⋅D/2 = ε|E|²/2, which is exactly the integral you want for the denominator, divided by 2.

[stevengj]

When comparing to a finite difference approximation dω/dR ≈ [ω(R+ΔR)–ω(R)] / ΔR, be aware that there are two convergence parameters: the resolution and ΔR.

In principle, to get agreement with perturbation theory, you need to take the limit as resolution→∞ and then take the limit as ΔR→0. One way to do this would be to keep doubling the resolution at a given ΔR until it is converged to your satisfaction, then halve ΔR and repeat, ….

[danielwboyce]

Been working some more on this. From the comments above I was able to get a calculation for dω/dR working, and then I starting comparing what perturbation theory would expect ω to be for a certain perturbation dR, and I got results that probably roughly match what we would expect (this is for the Ez case):

and we see error between what perturbation theory predicts for ω at R+dR and what Harminv finds for ω at R+dR increases following a power rule. In fact, it seems to follow a power rule too well and makes me think that I'm doing something wrong, but I haven't found anything concrete suggesting I've made a mistake. I also calculated the relative error between the value of dω/dR calculated with perturbation theory and what dω/dR on average would be (dω/dR = (ω(R+dR) - ω(R) / dR), but the strange thing is that error seems to decrease as dR increases.... Any ideas why that is?

My complete script is here

from __future__ import division

import meep as mp
import numpy as np
from statistics import mean
import matplotlib.pyplot as plt

def main():
    n = 3.4                 # index of waveguide
    r = 1
    a = r                   # inner radius of ring
    w = 1                   # width of waveguide
    b = a + w               # outer radius of ring
    pad = 4                 # padding between waveguide and edge of PML

    dpml = 2                # thickness of PML
    pml_layers = [mp.PML(dpml)]

    resolution = 100

    sr = b + pad + dpml            # radial size (cell is from 0 to sr)
    dimensions = mp.CYLINDRICAL    # coordinate system is (r,phi,z) instead of (x,y,z)
    cell = mp.Vector3(sr, 0, 0)

    m = 4

    geometry = [mp.Block(center=mp.Vector3(a + (w / 2)),
                         size=mp.Vector3(w, 1e20, 1e20),
                         material=mp.Medium(index=n))]

    # Finding a resonance mode with a high Q-value (calculated with Harminv)

    fcen = 0.15         # pulse center frequency
    df = 0.1            # pulse width (in frequency)

    sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df), mp.Ez, mp.Vector3(r+0.1))]

    sim = mp.Simulation(cell_size=cell,
                        geometry=geometry,
                        boundary_layers=pml_layers,
                        resolution=resolution,
                        sources=sources,
                        dimensions=dimensions,
                        m=m)

    h = mp.Harminv(mp.Ez, mp.Vector3(r+0.1), fcen, df)
    sim.run(mp.after_sources(h), until_after_sources=200)

    Harminv_freq_at_R = h.modes[0].freq

    sim.reset_meep()

    fcen = Harminv_freq_at_R
    df = 0.01

    sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df), mp.Ez, mp.Vector3(r + 0.1))]

    sim = mp.Simulation(cell_size=cell,
                        geometry=geometry,
                        boundary_layers=pml_layers,
                        resolution=resolution,
                        sources=sources,
                        dimensions=dimensions,
                        m=m)

    sim.run(until_after_sources=200)

    npts = 10
    angles = 2 * np.pi / npts * np.arange(npts)
    parallel_fields = []
    for angle in angles:
        point = mp.Vector3(a, angle)
        e_z_field = sim.get_field_point(mp.Ez, point)
        e_total_field = np.real(np.sqrt(e_z_field * np.conj(e_z_field)))
        parallel_fields.append(e_total_field)

        point = mp.Vector3(b, angle)
        e_z_field = sim.get_field_point(mp.Ez, point)
        e_total_field = np.real(np.sqrt(e_z_field * np.conj(e_z_field)))
        parallel_fields.append(e_total_field)

    numerator_surface_integral = 2 * np.pi * b * mean(parallel_fields)
    denominator_surface_integral = sim.electric_energy_in_box(center=mp.Vector3((b + pad/2) / 2), size=mp.Vector3(b + pad/2))
    perturb_theory_dw_dR = -Harminv_freq_at_R * numerator_surface_integral / (4 * denominator_surface_integral)

    center_diff_dw_dR = []
    Harminv_freqs_at_R_plus_dR = []

    drs = np.logspace(start=-3, stop=-1, num=10)

    for dr in drs:
        sim.reset_meep()
        w = 1 + dr  # width of waveguide
        b = a + w
        print(f'The current dr is dr={dr}')
        if len(Harminv_freqs_at_R_plus_dR) == 0:
            fcen = Harminv_freq_at_R
        else:
            fcen = Harminv_freqs_at_R_plus_dR[-1]
        df = 0.01

        sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df), mp.Ez, mp.Vector3(r + 0.1))]

        geometry = [mp.Block(center=mp.Vector3(a + (w / 2)),
                             size=mp.Vector3(w, 1e20, 1e20),
                             material=mp.Medium(index=n))]

        sim = mp.Simulation(cell_size=cell,
                            geometry=geometry,
                            boundary_layers=pml_layers,
                            resolution=resolution,
                            sources=sources,
                            dimensions=dimensions,
                            m=m)

        h = mp.Harminv(mp.Ez, mp.Vector3(r + 0.1), fcen, df)
        sim.run(mp.after_sources(h), until_after_sources=200)

        Harminv_freq_at_R_plus_dR = h.modes[0].freq
        Harminv_freqs_at_R_plus_dR.append(Harminv_freq_at_R_plus_dR)

        dw_dR = (Harminv_freq_at_R_plus_dR - Harminv_freq_at_R) / dr
        center_diff_dw_dR.append(dw_dR)

    relative_errors_dw_dR = [abs((dw_dR - perturb_theory_dw_dR) / dw_dR) for dw_dR in center_diff_dw_dR]

    perturb_predicted_freqs_at_R_plus_dR = [dr * perturb_theory_dw_dR + Harminv_freq_at_R for dr in drs]
    relative_errors_freqs_at_R_plus_dR = [abs((perturb_predicted_freqs_at_R_plus_dR[i] - Harminv_freqs_at_R_plus_dR[i]) / Harminv_freqs_at_R_plus_dR[i]) for i in range(len(Harminv_freqs_at_R_plus_dR))]

    if mp.am_master():
        plt.figure(dpi=150)
        plt.loglog(drs, relative_errors_dw_dR, 'bo-', label='relative error')
        plt.grid(True, which='both', ls='-')
        plt.xlabel('perturbation amount $dr$')
        plt.ylabel('relative error between $dω/dR$')
        plt.legend(loc='upper right')
        plt.title('Comparison of Perturbation Theory and \nCenter-Difference Calculations in Finding $dω/dR$')
        plt.tight_layout()
        # plt.show()
        plt.savefig('ring_Ez_perturbation_theory.dw_dR_error.png')
        plt.clf()

        plt.figure(dpi=150)
        plt.loglog(drs, relative_errors_freqs_at_R_plus_dR, 'bo-', label='relative error')
        plt.grid(True, which='both', ls='-')
        plt.xlabel('perturbation amount $dr$')
        plt.ylabel('relative error between $ω(R+dR)$')
        plt.legend(loc='upper left')
        plt.title('Comparison of resonance frequencies at $R+dR$ predicted by\nperturbation theory and found with Harminv')
        plt.tight_layout()
        # plt.show()
        plt.savefig('ring_Ez_perturbation_theory.freqs_error.png')
        plt.clf()

if __name__ == '__main__':
    main()

[danielwboyce]

My script for the Hz case was essentially the same, except in the part where I calculated the surface integrals the code was

    npts = 10
    angles = 2 * np.pi / npts * np.arange(npts)
    parallel_fields = []
    perpendicular_fields = []

    for angle in angles:
        point = mp.Vector3(a, angle)
        e_r_field = sim.get_field_point(mp.Er, point)
        temp_perpendicular_field = np.real(np.sqrt(e_r_field*np.conj(e_r_field)))
        perpendicular_fields.append(temp_perpendicular_field)

        e_p_field = sim.get_field_point(mp.Ep, point)
        e_z_field = sim.get_field_point(mp.Ez, point)
        temp_parallel_field = np.real(np.sqrt(e_p_field*np.conj(e_p_field) + e_z_field*np.conj(e_z_field)))
        parallel_fields.append(temp_parallel_field)

        point = mp.Vector3(b, angle)
        e_r_field = sim.get_field_point(mp.Er, point)
        temp_perpendicular_field = np.real(np.sqrt(e_r_field * np.conj(e_r_field)))
        perpendicular_fields.append(temp_perpendicular_field)

        e_p_field = sim.get_field_point(mp.Ep, point)
        e_z_field = sim.get_field_point(mp.Ez, point)
        temp_parallel_field = np.real(np.sqrt(e_p_field * np.conj(e_p_field) + e_z_field * np.conj(e_z_field)))
        parallel_fields.append(temp_parallel_field)

    numerator_surface_integral = 2 * np.pi * b * (mean(parallel_fields) - mean(perpendicular_fields))

which contrasts from the Ez case (to subtract the perpendicular components of the fields from the parallel ones). The error calculations for ω and dω/dR ended up looking almost the same as for the Ez case, as well: so if I made any mistakes, I at least made them consistently.

Any feedback that you may have on these results would be appreciated.

[oskooi]

It is more straightforward to separate the resonant mode into the two orthogonal polarizations (E_∥ and E_⟂) by setting up the simulation in Cartesian rather than Cylindrical coordinates. For a 2d cell in the xy plane, as long as k_z=0 (which is true in this example), depending on the current source Meep automatically separates the fields into two polarizations: (1) out-of-plane (H_x,H_y,E_z) and (2) in-plane (E_x,E_y,H_z).

[oskooi]

Also, the computation of the numerator (numerator_surface_integral) of the perturbation-theory formula is incorrect. Note e.g. that there is a missing (ε₁-ε₂) term which multiplies |E_∥|². For the inner ring surface, this term would be (1-n²) where n is the refractive index of the ring, while for the outer ring surface it is (n²-1). The key point is that you need to integrate correctly over all the perturbed fields in the structure.

[danielwboyce]

It is more straightforward to separate the resonant mode into the two orthogonal polarizations (E∥ and E⟂) by setting up the simulation in Cartesian rather than Cylindrical coordinates. For a 2d cell in the xy plane, as long as kz=0 (which is true in this example), depending on the current source Meep automatically separates the fields into two polarizations: (1) out-of-plane (Hx,Hy,Ez) and (2) in-plane (Ex,Ey,Hz).

Running the computation in cylindrical coordinates goes much master, though, especially at higher resolutions. Unless it's technically incorrect to use Ephi and Ez as the parallel E fields and Er as the perpendicular one, it would be much more efficient to continue to use cylindrical coordinates, though I'll edit my code to hopefully make it more readable (also to add the Δε and Δ(ε^-1) terms).

[stevengj]

use Ephi and Ez as the parallel E fields and Er as the perpendicular

This is correct, I agree. 👍

(As before, you can use an Ez source for the out-of-plane polarization, and an Hz source for the in-plane polarization.)

[stevengj]

Note that when computing the term, you should just use the D field. (Both ε and E⟂ are discontinuous, but D⟂ is continuous. You want Meep to compute D for you so that it handles the grid in the right way.)

[oskooi]

The following script, which is an extension/modification of the one from #1068:comment, computes dω/dR using perturbation theory and a finite(/forward) difference for either the E_∥ or E_⟂ fields.

The results for E_∥ are correct (i.e., the values produced by the two methods are comparable) but those for the perpendicular field are not (even when |εE_⟂|² is computed as |D_⟂|² where D_⟂=D_r). Since there does not seem to be any obvious bugs in the way the perturbation theory results are computed, is it possible that there might be a bug in electric_energy_in_box in cylindrical coordinates that is exposed only when using perpendicular fields?

import meep as mp
import numpy as np
import argparse

def main(args):
    if args.perpendicular:
        src_cmpt = mp.Hz
        fcen = 0.21        # pulse center frequency                                                                                                                                                                                                                                                   
    else:
        src_cmpt = mp.Ez
        fcen = 0.15        # pulse center frequency                                                                                                                           

    n = 3.4                 # index of waveguide                                                                                                                              
    w = 1                   # ring width                                                                                                                                      
    r = 1                   # inner radius of ring                                                                                                                            
    pad = 4                 # padding between waveguide and edge of PML                                                                                                       
    dpml = 2                # thickness of PML                                                                                                                                
    m = 5                   # angular dependence                                                                                                                              

    pml_layers = [mp.PML(dpml)]

    sr = r + w + pad + dpml        # radial size (cell is from 0 to sr)                                                                                                       
    dimensions = mp.CYLINDRICAL    # coordinate system is (r,phi,z) instead of (x,y,z)                                                                                        
    cell = mp.Vector3(sr)

    geometry = [mp.Block(center=mp.Vector3(r + (w / 2)),
                         size=mp.Vector3(w, mp.inf, mp.inf),
                         material=mp.Medium(index=n))]

    # find resonant mode with high Q using Harminv and broadband source                                                                                                       

    df = 0.1            # pulse width (in frequency)                                                                                                                          

    sources = [mp.Source(mp.GaussianSource(fcen,fwidth=df),
                         component=src_cmpt,
                         center=mp.Vector3(r+0.1))]

    sim = mp.Simulation(cell_size=cell,
                        geometry=geometry,
            boundary_layers=pml_layers,
                        resolution=args.res,
            sources=sources,
                        dimensions=dimensions,
            m=m)

    h = mp.Harminv(src_cmpt, mp.Vector3(r+0.1), fcen, df)
    sim.run(mp.after_sources(h),
            until_after_sources=100)

    frq_unperturbed = h.modes[0].freq

    sim.reset_meep()

    # unperturbed geometry with narrowband source centered at resonant frequency                                                                                              

    fcen = frq_unperturbed
    df = 0.05*fcen

    sources = [mp.Source(mp.GaussianSource(fcen,fwidth=df),
                         component=src_cmpt,
                         center=mp.Vector3(r+0.1))]

    sim = mp.Simulation(cell_size=cell,
                        geometry=geometry,
                        boundary_layers=pml_layers,
                        resolution=args.res,
                        sources=sources,
                        dimensions=dimensions,
                        m=m)

    sim.run(until_after_sources=100)

    if args.perpendicular:
        deps_inv = 1 - 1/n**2
        numerator_integral = deps_inv*2*np.pi*(-r*abs(sim.get_field_point(mp.Dr, mp.Vector3(r)))**2 + (r+w)*abs(sim.get_field_point(mp.Dr, mp.Vector3(r+w)))**2)
    else:
        deps = 1 - n**2
        numerator_integral = deps*2*np.pi*(r*abs(sim.get_field_point(mp.Ez, mp.Vector3(r)))**2 - (r+w)*abs(sim.get_field_point(mp.Ez, mp.Vector3(r+w)))**2)

    denominator_integral = sim.electric_energy_in_box(center=mp.Vector3(0.5*sr), size=mp.Vector3(sr))
    perturb_theory_dw_dR = -frq_unperturbed * numerator_integral / (4 * denominator_integral)

    # perturbed geometry with narrowband source                                                                                                                               

    dr = 0.01

    sim.reset_meep()

    sources = [mp.Source(mp.GaussianSource(fcen,fwidth=df),
                         component=src_cmpt,
                         center=mp.Vector3(r + dr + 0.1))]

    geometry = [mp.Block(center=mp.Vector3(r + dr + (w / 2)),
                         size=mp.Vector3(w, mp.inf, mp.inf),
                         material=mp.Medium(index=n))]

    sim = mp.Simulation(cell_size=cell,
                        geometry=geometry,
                        boundary_layers=pml_layers,
                        resolution=args.res,
                        sources=sources,
                        dimensions=dimensions,
                        m=m)

    h = mp.Harminv(src_cmpt, mp.Vector3(r+dr+0.1), fcen, df)
    sim.run(mp.after_sources(h),
            until_after_sources=100)

    frq_perturbed = h.modes[0].freq

    finite_diff_dw_dR = (frq_perturbed - frq_unperturbed) / dr

    print("dwdR:, {} (pert. theory), {} (finite diff.)".format(perturb_theory_dw_dR,finite_diff_dw_dR))

if __name__ == '__main__':
    parser = argparse.ArgumentParser()
    parser.add_argument('-perpendicular', action='store_true', help='use perpendicular field source (default: parallel field source)')
    parser.add_argument('-res', type=int, default=100, help='resolution (default: 100 pixels/um)')
    args = parser.parse_args()
    main(args)

Results for E_⟂ (via the -perpendicular runtime argument) at resolutions of 100 and 200 show the disparity:

resolution = 100

dwdR:, -0.012094433744995537 (pert. theory), -0.07980874733692356 (finite diff.)

resolution = 200

dwdR:, -0.011803783346264796 (pert. theory), -0.0798088157000526 (finite diff.)

The finite-difference result is larger than perturbation theory by a factor of 6.8.

Results for E_∥ at resolutions of 100 and 200, however, are comparable:

resolution = 100

dwdR:, -0.0854469622600482 (pert. theory), -0.08521213871394706 (finite diff.)

resolution = 100

dwdR:, -0.08544607152066104 (pert. theory), -0.085211185298259 (finite diff.)

For reference, the E_∥ mode has a frequency of 0.17 and Q of 1800 and the E_⟂ mode has 0.21 and 1200.

[danielwboyce]

It's necessary to calculate the phi-based E-field component in this calculation. It's treated as a parallel field.

import meep as mp
import numpy as np
import argparse

def main(args):
    n = 3.4                 # index of waveguide
    r = 1
    a = r                   # inner radius of ring
    w = 1                   # width of waveguide
    b = a + w               # outer radius of ring
    pad = 4                 # padding between waveguide and edge of PML

    dpml = 2                # thickness of PML
    pml_layers = [mp.PML(dpml)]

    resolution = args.res

    sr = b + pad + dpml            # radial size (cell is from 0 to sr)
    dimensions = mp.CYLINDRICAL    # coordinate system is (r,phi,z) instead of (x,y,z)
    cell = mp.Vector3(sr, 0, 0)

    geometry = [mp.Block(center=mp.Vector3(a + (w / 2)),
                         size=mp.Vector3(w, 1e20, 1e20),
                         material=mp.Medium(index=n))]

    # Finding a resonance mode with a high Q-value (calculated with Harminv)

    fcen = 0.18         # pulse center frequency
    df = 0.1            # pulse width (in frequency)

    m = 5

    if args.perp:
        component = mp.Hz
        component_name = 'meep.Hz'
    else:
        component = mp.Ez
        component_name = 'meep.Ez'

    sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df),
                         component=component,
                         center=mp.Vector3(r+0.1),
                         amplitude=1)]

    sim = mp.Simulation(cell_size=cell,
                        geometry=geometry,
                        boundary_layers=pml_layers,
                        resolution=resolution,
                        sources=sources,
                        dimensions=dimensions,
                        m=m)

    h = mp.Harminv(component, mp.Vector3(r+0.1), fcen, df)
    sim.run(mp.after_sources(h),
            until_after_sources=200)

    harminv_freq_at_r = h.modes[0].freq

    sim.reset_meep()

    # unperturbed geometry with narrowband source centered at resonant frequency

    fcen = harminv_freq_at_r
    df = 0.01

    sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df),
                         component=component,
                         center=mp.Vector3(r + 0.1),
                         amplitude=1)]

    sim = mp.Simulation(cell_size=cell,
                        geometry=geometry,
                        boundary_layers=pml_layers,
                        resolution=resolution,
                        sources=sources,
                        dimensions=dimensions,
                        m=m)

    sim.run(until_after_sources=200)

    # now need to calculate the surface integrals that go into dw/dR. Fields parallel and perpendicular to the interface
    # AND at the inner and outer surfaces are treated differently, so each will be calculated separately.

    angle = 0
    deps_inner = 1 - n ** 2
    deps_inv_inner = 1 - 1/(n**2)
    deps_outer = n ** 2 - 1
    deps_inv_outer = -1 + 1 / (n ** 2)

    # section for fields at inner surface
    point = mp.Vector3(a, angle)

    # section for fields parallel to interface (Ez and Ep)
    e_z_field = abs(sim.get_field_point(mp.Ez, point))**2
    e_p_field = abs(sim.get_field_point(mp.Ep, point))**2
    parallel_field = e_z_field + e_p_field
    # fields have to be multiplied by Δε to get integrand
    parallel_field_integrand = deps_inner * parallel_field
    parallel_fields_inner = parallel_field_integrand * 2 * np.pi * a

    # section for fields perpendicular to interface (Dr)
    d_r_field = abs(sim.get_field_point(mp.Dr, point))**2
    perpendicular_field = d_r_field
    # fields have to be multiplied by Δ(1/ε) to get integrand
    perpendicular_field_integrand = deps_inv_inner * perpendicular_field
    perpendicular_fields_inner = perpendicular_field_integrand * 2 * np.pi * a

    # section for fields at outer surface
    point = mp.Vector3(b, angle)

    # section for fields parallel to interface (Ez and Ep)
    e_z_field = abs(sim.get_field_point(mp.Ez, point))**2
    e_p_field = abs(sim.get_field_point(mp.Ep, point))**2
    parallel_field = e_z_field + e_p_field
    # fields have to be multiplied by Δε to get integrand
    parallel_field_integrand = deps_outer * parallel_field
    parallel_fields_outer = parallel_field_integrand * 2 * np.pi * b

    # section for fields perpendicular to interface (Dr)
    d_r_field = abs(sim.get_field_point(mp.Dr, point))
    perpendicular_field = d_r_field**2
    # fields have to be multiplied by Δ(1/ε) to get integrand
    perpendicular_field_integrand = deps_inv_outer * perpendicular_field
    perpendicular_fields_outer = perpendicular_field_integrand * 2 * np.pi * b

    numerator_surface_integral = parallel_fields_outer - perpendicular_fields_outer + parallel_fields_inner - perpendicular_fields_inner
    denominator_surface_integral = sim.electric_energy_in_box(center=mp.Vector3(), size=mp.Vector3(2*sr))
    perturb_theory_dw_dr = -harminv_freq_at_r * numerator_surface_integral / (2 * denominator_surface_integral)

    # perturbed geometry with narrowband source

    dr = args.dr

    sim.reset_meep()
    a = r + dr  # inner radius of ring
    b = a + w  # outer radius of ring

    fcen = perturb_theory_dw_dr * dr + harminv_freq_at_r
    df = 0.01

    sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df),
                         component=component,
                         center=mp.Vector3(r + dr + 0.1))]

    geometry = [mp.Block(center=mp.Vector3(a + (w / 2)),
                         size=mp.Vector3(w, 1e20, 1e20),
                         material=mp.Medium(index=n))]

    sim = mp.Simulation(cell_size=cell,
                        geometry=geometry,
                        boundary_layers=pml_layers,
                        resolution=resolution,
                        sources=sources,
                        dimensions=dimensions,
                        m=m)

    h = mp.Harminv(component, mp.Vector3(r + dr + 0.1), fcen, df)
    sim.run(mp.after_sources(h),
            until_after_sources=200)

    harminv_freq_at_r_plus_dr = h.modes[0].freq

    finite_diff_dw_dr = (harminv_freq_at_r_plus_dr - harminv_freq_at_r) / dr

    print("component:, {}".format(component_name))
    print("res:, {}".format(resolution))
    print("dr:, {}".format(dr))
    print("w:, {} (unperturbed), {} (perturbed)".format(harminv_freq_at_r, harminv_freq_at_r_plus_dr))
    print("dwdR:, {} (pert. theory), {} (finite diff.)".format(perturb_theory_dw_dr, finite_diff_dw_dr))

if __name__ == '__main__':
    parser = argparse.ArgumentParser()
    parser.add_argument('-perp', type=bool, default=False, help='controls if in-plane perpendicular fields are excited or not. (default: False, for Ez source)')
    parser.add_argument('-res', type=int, default=100, help='resolution (default: 100 pixels/um)')
    parser.add_argument('-dr', type=float, default=0.02, help='dr (default: 0.02)')
    args = parser.parse_args()
    main(args)

For $dr=0.01$, the results at resolutions 100 and 200 in the perpendicular field case show agreement within 1%

resolution = 100

dwdR:, -0.08050386030223725 (pert. theory), -0.07980868621023374 (finite diff.)

resolution = 200

dwdR:, -0.08020283538083425 (pert. theory), -0.07980875403289789 (finite diff.)

Results for the parallel case are identical to before

resolution = 100

dwdR:, -0.08544696870452953 (pert. theory), -0.08521245606197825 (finite diff.)

resolution = 200

dwdR:, -0.08544608457833698 (pert. theory), -0.08521149913008619 (finite diff.)

[oskooi]

Yes, you are correct. I was just about to comment on this with the same results that you posted. As you pointed out, the correct way to compute the numerator for the H_z source must include D_r and E_φ:

    deps = 1 - n**2
    deps_inv = 1 - 1/n**2

    if args.perpendicular:
    para_integral = deps*2*np.pi*(r*abs(sim.get_field_point(mp.Ep, mp.Vector3(r)))**2 - (r+w)*abs(sim.get_field_point(mp.Ep, mp.Vector3(r+w)))**2)
        perp_integral = deps_inv*2*np.pi*(-r*abs(sim.get_field_point(mp.Dr, mp.Vector3(r)))**2 + (r+w)*abs(sim.get_field_point(mp.Dr, mp.Vector3(r+w)))**2)
        numerator_integral = para_integral + perp_integral
    else:
        numerator_integral = deps*2*np.pi*(r*abs(sim.get_field_point(mp.Ez, mp.Vector3(r)))**2 - (r+w)*abs(sim.get_field_point(mp.Ez, mp.Vector3(r+w)))**2)

[oskooi]

Also, for reference, the same calculation in 2d Cartesian coordinates which we may or may not want to include in the tutorial only to demonstrate the performance disparity with cylindrical coordinates:

import meep as mp
import numpy as np
import argparse

def main(args):
    if args.perpendicular:
    src_cmpt = mp.Hz
    fcen = 0.21         # pulse center frequency                                                                                                                              
    else:
        src_cmpt = mp.Ez
        fcen = 0.17         # pulse center frequency                                                                                                                              

    n = 3.4                 # index of waveguide                                                                                                                                  
    w = 1                   # ring width                                                                                                                                          
    r = 1                   # inner radius of ring                                                                                                                                
    pad = 4                 # padding between waveguide and edge of PML                                                                                                           
    dpml = 2                # thickness of PML                                                                                                                                    

    pml_layers = [mp.PML(dpml)]

    sxy = 2*(r+w+pad+dpml)
    cell_size = mp.Vector3(sxy,sxy)

    symmetries = [mp.Mirror(mp.X,phase=+1 if args.perpendicular else -1),
                  mp.Mirror(mp.Y,phase=-1 if args.perpendicular else +1)]

    geometry = [mp.Cylinder(material=mp.Medium(index=n),
                            radius=r+w,
                            height=mp.inf),
                mp.Cylinder(material=mp.vacuum,
                            radius=r,
                            height=mp.inf)]

    # find resonant mode with high Q using Harminv and broadband source                                                                                                           

    df = 0.2*fcen            # pulse width (in frequency)                                                                                                                         

    sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df),
                         component=src_cmpt,
                         center=mp.Vector3(r+0.1)),
           mp.Source(mp.GaussianSource(fcen, fwidth=df),
                         component=src_cmpt,
                         center=mp.Vector3(-(r+0.1)),
                         amplitude=-1)]

    sim = mp.Simulation(cell_size=cell_size,
            geometry=geometry,
            boundary_layers=pml_layers,
            resolution=args.res,
            sources=sources,
            symmetries=symmetries)

    h = mp.Harminv(src_cmpt, mp.Vector3(r+0.1), fcen, df)
    sim.run(mp.after_sources(h),
            until_after_sources=100)

    frq_unperturbed = h.modes[0].freq

    sim.reset_meep()

    # unperturbed geometry with narrowband source centered at resonant frequency                                                                                                  

    fcen = frq_unperturbed
    df = 0.05*fcen

    sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df),
                         component=src_cmpt,
                         center=mp.Vector3(r+0.1)),
               mp.Source(mp.GaussianSource(fcen, fwidth=df),
                         component=src_cmpt,
                         center=mp.Vector3(-(r+0.1)),
                         amplitude=-1)]

    sim = mp.Simulation(cell_size=cell_size,
                        geometry=geometry,
                        boundary_layers=pml_layers,
                        resolution=args.res,
                        sources=sources,
                        symmetries=symmetries)

    sim.run(until_after_sources=100)

    deps = 1 - n**2
    deps_inv = 1 - 1/n**2

    if args.perpendicular:
        para_integral = deps*2*np.pi*(r*abs(sim.get_field_point(mp.Ey, mp.Vector3(r)))**2 - (r+w)*abs(sim.get_field_point(mp.Ey, mp.Vector3(r+w)))**2)
        perp_integral = deps_inv*2*np.pi*(-r*abs(sim.get_field_point(mp.Dy, mp.Vector3(y=r)))**2 + (r+w)*abs(sim.get_field_point(mp.Dy, mp.Vector3(y=r+w)))**2)
    numerator_integral = para_integral + perp_integral
    else:
        numerator_integral = deps*2*np.pi*(r*abs(sim.get_field_point(mp.Ez, mp.Vector3(r)))**2 - (r+w)*abs(sim.get_field_point(mp.Ez, mp.Vector3(r+w)))**2)

    denominator_integral = sim.electric_energy_in_box(center=mp.Vector3(), size=mp.Vector3(sxy-2*dpml,sxy-2*dpml))
    perturb_theory_dw_dR = -frq_unperturbed * numerator_integral / (8 * denominator_integral)

    # perturbed geometry with narrowband source                                                                                                                                   

    dr = 0.04

    sim.reset_meep()

    sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df),
                         component=src_cmpt,
                         center=mp.Vector3(r+dr+0.1)),
               mp.Source(mp.GaussianSource(fcen, fwidth=df),
                         component=src_cmpt,
                         center=mp.Vector3(-(r+dr+0.1)),
                         amplitude=-1)]

    geometry = [mp.Cylinder(material=mp.Medium(index=n),
                            radius=r+dr+w,
                            height=mp.inf),
                mp.Cylinder(material=mp.vacuum,
                            radius=r+dr,
                            height=mp.inf)]

    sim = mp.Simulation(cell_size=cell_size,
                        geometry=geometry,
                        boundary_layers=pml_layers,
                        resolution=args.res,
                        sources=sources,
                        symmetries=symmetries)

    h = mp.Harminv(src_cmpt, mp.Vector3(r+dr+0.1), fcen, df)
    sim.run(mp.after_sources(h),
            until_after_sources=100)

    frq_perturbed = h.modes[0].freq

    finite_diff_dw_dR = (frq_perturbed - frq_unperturbed) / dr

    print("dwdR:, {} (pert. theory), {} (finite diff.)".format(perturb_theory_dw_dR,finite_diff_dw_dR))

if __name__ == '__main__':
    parser = argparse.ArgumentParser()
    parser.add_argument('-perpendicular', action='store_true', help='use perpendicular field source (default: parallel field source)')
    parser.add_argument('-res', type=int, default=30, help='resolution (default: 30 pixels/um)')
    args = parser.parse_args()
    main(args)

The two results are again comparable and similar to those for cylindrical coordinates but the resolution is much lower (30 vs. 100/200). This is because the 2d computation takes longer even with the two mirror-symmetry planes which reduce the cell size by a factor of 4.

H_z source:

resolution = 30

dwdR:, -0.08730321280915465 (pert. theory), -0.07964697654790631 (finite diff.)

E_z source:

resolution = 30

dwdR:, -0.08537783093554978 (pert. theory), -0.08379133012667156 (finite diff.)

[stevengj]

Note that in the Cartesian-coordinate (2d) implementation, the mode produced by a point source is actually the cos(mφ) mode, not exp(imφ), or equivalently it is the superposition of the ±m modes. This means that when you do the numerator integral, it is not simply 2πr times the value of the point — it is the integral of cos²(mφ), which gives a factor of 1/2. (Of course, for general non-circular shapes in 2d you would have to do the surface integral numerically.)