Tutorial for nonaxisymmetric point-dipole sources in cylindrical coordinates

[oskooi]

Closes #2108.

This is an initial and incomplete draft intended for early review and feedback.

[smartalecH]

This looks great, @oskooi. Nice job!

The idea here, of course, is that we can efficiently simulate the effects of a single dipole source in 3D using a few (or more than a few) relatively cheap cylindrical simulations. As we've discussed, it may not actually be that much cheaper (but hopefully it is, especially since it parallelizes much better!)

It might be worth it to have a bit of dialogue discussing this nuance. And maybe link to the other tutorial which exploits linear time invariance to compute the LEE of an entire structure using several dipoles. In other words, you can do various simulations at multiple r positions (each with multiple m numbers) and compute the LEE for the whole 3D device.

That naturally begs the question: how easily we can adapt the other LEE tutorial (using custom basis functions) within cylindrical coordinates?

Finally, (as another todo) we could add a tutorial that does relatively simple adjoint optimization on the LEE. These results, together with @mochen4's LDOS adjoint should enable that (in fact, the LDOS saves us from having to do an adjoint solve, right, Mo?) There would be a lot going on with this kind of thing... but it would be rather interesting!

[oskooi]

I modified the tutorial to include results which validate the scaling factors necessary to convert the radiated flux from the LED simulation in cylindrical coordinates from a point dipole at $r = 0$ and $r > 0$ to 3d. The agreement between cylindrical coordinates and 3d is good. (We may want to consider adapting this tutorial into a unit test.) What is still missing is a first-principles explanation for these scaling factors.

As described in #2108 (comment), it seems computing the flux for an $r > 0$ point dipole using a Fourier-series expansion is only accurate when the point source is sufficiently far away from $r = 0$. It's not clear whether this is a feature or a bug in the way we have set up the computation.

[mochen4]

Finally, (as another todo) we could add a tutorial that does relatively simple adjoint optimization on the LEE. These results, together with @mochen4's LDOS adjoint should enable that (in fact, the LDOS saves us from having to do an adjoint solve, right, Mo?) There would be a lot going on with this kind of thing... but it would be rather interesting!

Actually, we still need an adjoint solve for the LDOS, because the objective can be a function of LDOS with several sources at different frequencies. We need to take the derivatives for adjoint sources and adjoint runs. If the objective is just one LDOS itself, then in principle we could save the adjoint solve.

[stevengj]

I don't like putting an empirical scale factor in the tutorial. If we want to claim this is correct, someone should go through the algebra and derive it — that pi^5 factor looks weird. The scale factor for r=0 is much simpler to explain (though even then more work is required to get the factor of 2).

Better to keep the tutorial focused on dimensionless quantities like extraction efficiency or Purcell enhancement so that these scale factors are irrelevant.

[smartalecH]

We mentioned that one possible reason for the discrepancy w.r.t source position (along r) could be related to how we are placing the (delta) source.

A few thoughts to add to this: we are currently placing the source in cylindrical coordinates the same as we do in cartesian. That is, we use a restriction, which is the transpose of a linear interpolation. However, this same restriction (i.e. using the same code) in cylindrical coordinates will place an unequal weight for different points along r (the z direction restriction should be fine... it's just the r direction).

We may need to scale the restriction using something like the bilinear transform (or maybe just transform the restriction using the Jacobian of the coordinate transform?).

[stevengj]

@smartalecH, yes — see

meep/src/loop_in_chunks.cpp

Lines 112 to 113 in 1fe3899

	Integration Weights in Cylindrical Coordinates
	FIXME: implement this below?

[codecov-commenter]

Codecov Report

Merging #2452 (fdca861) into master (75ffb03) will increase coverage by 0.24%.
The diff coverage is n/a.

📣 This organization is not using Codecov’s GitHub App Integration. We recommend you install it so Codecov can continue to function properly for your repositories. Learn more

@@            Coverage Diff             @@
##           master    #2452      +/-   ##
==========================================
+ Coverage   72.66%   72.90%   +0.24%     
==========================================
  Files          17       17              
  Lines        5162     5212      +50     
==========================================
+ Hits         3751     3800      +49     
- Misses       1411     1412       +1     

see 1 file with indirect coverage changes

[oskooi]

Better to keep the tutorial focused on dimensionless quantities like extraction efficiency or Purcell enhancement so that these scale factors are irrelevant.

I have replaced the calculation of the dimensionful radiated flux with the dimensionless extraction efficiency.

Another useful demonstration that we may want to consider adding to this tutorial is computing the extraction efficiency for some narrow cone of angles (rather than all angles). See related discussion in #2368.

[smartalecH]

Another useful demonstration that we may want to consider adding to this tutorial is computing the extraction efficiency for some narrow cone of angles (rather than all angles).

Yes! As mentioned in that discussion, a classic near2far makes this rather trivial. However, Our current near2far is overly expensive for this kind of calculation. The current implementation uses the Green's function to project to any arbitrary point in space. But what we care about (in this particular case) is the result as function of angle.

Now you could use some trig to convert your near2far response so long as it's sufficiently far... but this is needlessly complicated. When evaluating infinitely far away, the Green's function simplifies significantly resulting in a simple Fourier transform integral (with some nuances in cylindrical coordinates...). This integral is not only cheaper to evaluate, but directly provides the far field as a function on angle (no additional post processing). You could trivially implement this in python even (as done here).

[smartalecH]

We mentioned that one possible reason for the discrepancy w.r.t source position (along r) could be related to how we are placing the (delta) source.

Here's an idea: rather than place a single dipole source, why not place a line source along the waveguide (really a circular plane source)? This way, only one r weight is needed, and it's far from the center anyway (so its error will be minimal).

Since we are already computing extraction efficiency, we might as well demonstrate it with the basis function approach (method 3) in the tutorials. For this case though, just simulate a single basis function (the DC term). The computational cost is the same as the current tutorial (for both the cylindrical and 3D case) but now the error is effectively eliminated.

If we have a discrepancy between the cylindrical and 3D case, then we know the primary error is not with the interpolation weights, but something else (although those weights should still be fixed eventually)

[stevengj]

Note that for an r=0 source, I don't think it should matter (modulo factors of 2) whether you do an rhat source or an rhat+iφhat source.

[oskooi]

Note that for an r=0 source, I don't think it should matter (modulo factors of 2) whether you do an rhat source or an rhat+iφhat source.

Using the necessary bug fix in #2459, I verified for the case of a point source at $r=0$ with $m=1$ that the two polarizations $\hat{r}$ and $\hat{r}+i\hat{\phi}$ produce the same result for the extraction efficiency.