Closed stevengj closed 4 years ago
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;
}
In particular, change get_eigenmode
to:
call ishomogeneous
when the maxwell data is initialized to check whether it is a homogeneous region.
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.)
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.
Note, by the way, that our k units already include the 2π factor, so you should use k = sqrt(ω²n² - (kₓ+m/Λ)²)
This should be fast enough now and gives zero for the evanescent modes.
We should probably re-open this because of
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.
Effectively, this is the inverse of the get_dominant_planewave
function (https://github.com/stevengj/mpb/pull/71): a set_dominant_planewave
function
The MPB function maxwell_set_planewave
was added by https://github.com/stevengj/mpb/commit/4b8d7a126ff0203fe627ada406793c59f1a89720 … now we just need a Meep interface to allow you to use this from the mode solver instead of supplying a band index.
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:
g[3]
input above. (Note that it will only need to compute 1 band in this case regardless of the order.)kperp = sqrt(n²ω² - (kparallel + G)²)
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.
@smartalecH, what API does Lumerical (or other commercial codes) expose for specifying planewave diffraction orders and amplitudes/polarizations?
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.
Is there documentation online you can link? I'm wondering how precisely you specify the polarization.
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.
For now we can just pass k
, a triplet g
of integers specifying the G
vector (diffraction order), an axis
vector, and the s
and p
amplitudes (or a flag indicating s or p) — these are the inputs to MPB's function. (s and p polarizations are defined relative to the (k+G)×axis
plane.)
Later on we can worry about a default axis and other niceties.
(In 2d, the axis
vector should probably default to z
?)
Things to do:
Modify fields::get_eigenmode
to accept optional additional parameters int g[3], double axis[3], std::complex<double> s, std::complex<double> p
. If band_num == 0
, then it will pass these to maxwell_set_planewave
as above instead of calling the eigensolver. Similarly for functions that call get_eigenmode
. (Alternatively, we might define a diffractedplanewave
class to encapsulate these parameters in a single object, and a pointer to this object, which defaults to NULL
.)
if (0)
in my patch above would instead be if (diffraction argument is present)
. And the axis
, g
, s
, and p
parameters that are hard-coded in the patch above would instead be taken from the diffractedplanewave *dp
object)Define something like a DiffractedPlanewave(g=(i,j,k), axis=(x,y,z), s=..., p=....)
class in Python to encapsulate these parameters.
Modify the Python eigenmode source and coefficients to optionally take a DiffractedPlanewave
object instead of a band number. If such an object is passed, then we pass it along to the C++ API above.
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