Optimize mode coefficients for planewave case

[stevengj]

For the mode-coefficient expansion, there is a special case worth optimizing. In particular, if:

• The flux plane is in a homogeneous medium. (Initially, I would just support homogeneous scalar ε and μ, with no anisotropy.)

• The flux plane fills the entire unit cell in those directions, and the boundary conditions in those directions are Bloch-periodic.

In this case, the propagating modes are a finite set of planewaves. Computing the mode coefficients is just a discrete Fourier transform. So one option would be to circumvent MPB entirely in this case and just perform a DFT to get the coefficients of the planewaves, appropriately normalized. It could even compute the coefficients of the evanescent waves too.

Alternatively, since MPB is also super-efficient in this case there is no huge need to optimize the eigenmode calculation. What would be nice, however, would be to:

• Automatically compute the correct number of (non-evanescent) modes for a given frequency.

• Give MPB an initial "guess" for each mode consisting of the analytically known k.

For example suppose we have a 2d calculation with period Λ in the x direction, bloch-periodic boundary conditions with a given kₓ, and have a flux plane of width Λ in a homogeneous region of index n. Then the propagating modes at a frequency ω (in c=1 units) are the real solutions of k = sqrt(ω²n² - (kₓ+2πm/Λ)²) for integers m. That is, you just (a) compute the integers m for which the square root is real = number of modes; and (b) compute the corresponding k values, passing (kₓ,k) to MPB to get the mode at that frequency.

cc @HomerReid

[stevengj]

In mpb.cpp, after you call set_maxwell_dielectric, you can check for a uniform material by something like:

bool ishomogeneous(maxwell_data *mdata) {
    const symmetric_matrix *eps_inv = mdata->eps_inv;
    int n = mdata->fft_output_size;
    for (int i = 0; i < n; ++i) if (eps_inv[i] != eps_inv[0]) return false;
    return true;
}

[stevengj]

In particular, change get_eigenmode to:

1. call ishomogeneous when the maxwell data is initialized to check whether it is a homogeneous region.

2. Return NULL if the requested band_num can't be found. This happens when match_frequency is true and either the Newton iteration fails to converge (which currently aborts) or it is homogeneous and band_num exceeds the analytically known number of planewave modes. The caller can then decide what to do — eigenmode sources will abort, but getting the eigenmode coefficient will simply return 0. (If the mode doesn't exist, there are still evanescent modes, but they carry no power.)

3. If it is homogeneous and match_frequency is true, then use the analytical formula to set the initial value of kpoint in the direction d. (Remember, this component of kpoint gets overwritten anyway by the Newton iteration in that case, so this is just improving the initial guess.)

This way the homogeneous case has exactly the same semantics as the inhomogeneous case, it is just faster and more reliable.

[stevengj]

Note, by the way, that our k units already include the 2π factor, so you should use k = sqrt(ω²n² - (kₓ+m/Λ)²)

[stevengj]

This should be fast enough now and gives zero for the evanescent modes.

[stevengj]

We should probably re-open this because of
#604: we'd like to ensure that it gets power-orthogonal modes in the 3d case and picks the polarization in a sensible deterministic way.

The easiest option here might be just to set the H matrix (which = modes in a transverse planewave basis) to have one coefficient = 1 and the other coefficients = 0, which gives a planewave mode. Then we can order the modes however we want. A little thought may be required to ensure that the modes are power-orthogonal, which we can do by ensuring that we choose modes in the s and p polarizations with respect to the incident plane.

[stevengj]

Effectively, this is the inverse of the get_dominant_planewave function (NanoComp/mpb#71): a set_dominant_planewave function

[stevengj]

The MPB function maxwell_set_planewave was added by NanoComp/mpb@4b8d7a1 … now we just need a Meep interface to allow you to use this from the mode solver instead of supplying a band index.

[stevengj]

Okay, set_planewave seems to be working. I tried it with the examples/oblique-planewave.py example, and simply hacked mpb.cpp to call set_planewave instead of solving the eigenproblem:

--- a/src/mpb.cpp
+++ b/src/mpb.cpp
@@ -434,7 +434,7 @@ void *fields::get_eigenmode(double omega_src, direction d, const volume where, c
   /*- part 2: newton iteration loop with call to MPB on each step */
   /*-         until eigenmode converged to requested tolerance    */
   /*--------------------------------------------------------------*/
-  if (am_master()) do {
+  if (0) do {
       eigensolver(H, eigvals, maxwell_operator, (void *)mdata,
 #if MPB_VERSION_MAJOR > 1 || (MPB_VERSION_MAJOR == 1 && MPB_VERSION_MINOR >= 6)
                   NULL, NULL, /* eventually, we can support mu here */
@@ -485,6 +485,30 @@ void *fields::get_eigenmode(double omega_src, direction d, const volume where, c
     } while (match_frequency &&
              fabs(sqrt(eigvals[band_num - 1]) - omega_src) > omega_src * match_tol);
 
+  {
+    int g[3] = {0,0,0};
+    mpb_real axis[3] = {0,1,0};
+    scalar_complex s = {1.0,0.0};
+    scalar_complex p = {0.0,0.0};
+
+    master_printf("DEBUG k = %g, %g, %g for omega = %g\n", k[0],k[1],k[2], omega_src);
+    update_maxwell_data_k(mdata, k, G[0], G[1], G[2]);
+    maxwell_set_planewave(mdata, H, band_num, g, s, p, axis);
+
+    eigvals[band_num - 1] = omega_src;
+
+    evectmatrix_resize(&W[0], 1, 0);
+    evectmatrix_resize(&W[1], 1, 0);
+    for (int i = 0; i < H.n; ++i)
+      W[0].data[i] = H.data[H.p - 1 + i * H.p];
+    maxwell_ucross_op(W[0], W[1], mdata, kdir); // W[1] = (dTheta/dk) W[0]
+    mpb_real vscratch;
+    evectmatrix_XtY_diag_real(W[0], W[1], &vgrp, &vscratch);
+    vgrp /= sqrt(eigvals[band_num - 1]);
+
+    master_printf("DEBUG got vgrp = %g\n", vgrp);
+  }
+
   double eigval = eigvals[band_num - 1];

and it produces the correct output.

Now the question is what API do we expose. Basically, instead of giving the band number, the user would specify:

1. The diffraction order in each direction (transverse to the source plane). This is the g[3] input above. (Note that it will only need to compute 1 band in this case regardless of the order.)
2. The frequency. This would be used to solve analytically for the k component normal to the source plane by the usual formula kperp = sqrt(n²ω² - (kparallel + G)²)
3. The complex amplitudes of the s and p polarizations relative to the incident plane, which is the plane orthogonal to the source plane and containing the k vector.

The k vector transverse to the source plane would be determined as it is now (from the boundary conditions etc). The only question is how to define the incident plane when k=0 (for the purpose of defining the s and p polarizations). We could choose it arbitrarily, of course, but it would be nice for the user to have a way to specify it so that they don't get a discontinuity when they change k from 0 to nonzero.

[stevengj]

@smartalecH, what API does Lumerical (or other commercial codes) expose for specifying planewave diffraction orders and amplitudes/polarizations?

[smartalecH]

To examine diffraction orders, Lumerical's FDTD uses a transmission_grating object, which is a modified DFTMonitor. It has some cool features, like spectral averaging (e.g. Welch's method), temporal apodization, and spatial interpolation.

There are various planewave and mode source objects that allow the user to specify the angle, amplitude, polarization/mode, etc. To launch a particular diffraction order, for example, you can specify a planewave with the corresponding angle, amplitude profile, and polarization using the planewave object.

[stevengj]

Is there documentation online you can link? I'm wondering how precisely you specify the polarization.

[smartalecH]

Here is some information about the plane source.

Here is some information about the DFT monitor used to extract the diffraction orders.

The online documentation is a little bare. I usually learn new features through various application examples. Here's a transmission grating example that uses the features you were asking about.

The source object only asks the user to specify the polarization angle, but not the actual polarization mode. The DFT monitors don't appear to do any mode decomposition on the Fourier transformed fields, so they don't care about the polarization either.