tutorial for oblique planewave in cylindrical coordinate

[mochen4]

Documentation and tutorial for how to have oblique incident planewave in cylindrical coordinate. Compared the scattering flux of a sphere under different angles of incidence. Closes #2630.

The setup of the tutorial looks like

Part of the source around $r=0$ is currently excluded for $|m|>1$ due to #2644. In addition, as noted in #2538 (comment), source at $r=0$ for $|m|=1$ may also not be accurate. On the other hand, both issues should only introduce 1st order errors, which should go away with resolution.

Here is a table of the results at various incidence angles at various cutoff of |m| at resolution=100.

Angle of incidence	$$|m|=0$$	$$|m| \le 1$$	$$|m| \le 2$$	$$|m| \le 3$$	$$|m| \le 4$$
$$0^\circ$$	0.	17.83312867	17.83312867	17.83312867	17.83312867
$$5^\circ$$	0.26219472	16.28351528	16.8657832	16.99266891	17.02334217
$$7.5^\circ$$	0.32401024	15.14771059	16.08409591	16.12975925	16.17413882
$$15^\circ$$	1.16801433	12.15571072	14.8310367	15.03711602	15.13986574

A higher cutoff should have better convergence, but as mentioned in #2644 (comment), a higher cutoff currently required excluding more pixels around $r=0$.

[oskooi]

Thanks for putting this together.

It would be good to provide some comments on the rate of convergence of the result with $m$ in the series expansion for increasing angles of the incident planewave. We could cite the analysis from a different tutorial.

[oskooi]

The comment should reference the bug reports #2538 and #2644 for the $E_r$ source at $r = 0$ for |m| > 0.

[mochen4]

I wonder whether we should try to fix #2538 or #2644 first, since It seems that they really affect the accuracy. I was hoping maybe we could fix some bugs and I can update the tutorial later; that's why I didn't reference the issues.

[oskooi]

Based on the workaround described in #2704, would you try using src_offset = 1.5 / resolution? This should hopefully improve the accuracy of your results.

[mochen4]

I tried that before but there wasn't much improvement. Due to #2644, additional pixels have to be excluded for |m|>1, and the larger |m| is, the more pixels needed to be excluded, which also implies that having more terms in the expansion may not help.

[mochen4]

Ensuring stability(#2644) and forcing complex fields (inspired by #2726 (comment)) produce much better agreement:

Angle of incidence	$$|m|=0$$	$$|m| \le 1$$	$$|m| \le 2$$	$$|m| \le 3$$	$$|m| \le 4$$
$$0^\circ$$	0.	8.89879305	8.89879305	8.89879305	8.89879305
$$5^\circ$$	0.38449451	8.43871285	8.71133751	8.71249049	8.7718544
$$7.5^\circ$$	0.66191608	8.28861876	8.71832748	8.72696144	8.75597351
$$15^\circ$$	1.84922004	7.4380022	8.64164363	8.75756958	8.79085675

Also notice that the size of source is increased to more accurately represent a planewave, since we don't have real planewave in cylindrical coordinate:

[oskooi]

Good to see this is finally working as expected. A few comments and questions:

• don't the angles from alpha_list = [0, np.pi/36, np.pi/24, np.pi/12] in the script correspond to 0°, 5°, 7.5°, and 15° rather than 0°, 5°, 10°, 15° in your table?
• why choose a maximum value of m of M = 4 for each angle? should M be a constant independent of the angle of the planewave?
• in the Jacobi-Anger expansion, is it necessary to perform a sum of m from -M to M? why not sum from 0 to M where the results for the flux for the m > 0 runs are multiplied by two? (Running the script on my machine with 2 Intel Xeon cores at 4.2 GHz takes ~40 minutes.)
• it would be useful to explain why it is necessary to set Courant = 0.2 for all the simulations in the expansion. Does the Courant factor need to be a constant? What about making its value a function of m?
• should the line source extend into the PML to the edge of the cell?

[mochen4]

Good to see this is finally working as expected. A few comments and questions:

• don't the angles from alpha_list = [0, np.pi/36, np.pi/24, np.pi/12] in the script correspond to 0°, 5°, 7.5°, and 15° rather than 0°, 5°, 10°, 15° in your table?

Thanks! I fixed the value.

• why choose a maximum value of m of M = 4 for each angle? should M be a constant independent of the angle of the planewave?

I chose M=4 for simplicity. I don't think larger angles would require more terms to converge, but I am not sure.

• in the Jacobi-Anger expansion, is it necessary to perform a sum of m from -M to M? why not sum from 0 to M where the results for the flux for the m > 0 runs are multiplied by two? (Running the script on my machine with 2 Intel Xeon cores at 4.2 GHz takes ~40 minutes.)

There might be a correspondence between m and -m terms; I will go through the math and make sure I have the correct signs and factor.

• it would be useful to explain why it is necessary to set Courant = 0.2 for all the simulations in the expansion. Does the Courant factor need to be a constant? What about making its value a function of m?

Thanks! I can make it as a function of m.

• should the line source extend into the PML to the edge of the cell?

For Cartesian cells, sources extending into the PML with the periodic BC will generate planewave; but in cylindrical coordinate that is not the case. I am not sure which way would be better.

[stevengj]

In the radial direction, it would be good to make the computation converge faster with the cell radius R so that you can use a smaller cell. A basic issue is that you are just truncating the current source at R, so you are introducing an error that converges as 1/R.

One way to get faster convergence would be to make the source decay smoothly in R, similar to the windowed-Green's function approach (#1952), e.g. by multiplying it by a bump function.

(The trick of extending the source into the PML and using periodic boundary conditions, which we do in Cartesian coordinates, doesn't seem applicable here since the problem cannot be periodic in r. Although maybe for the planewave you could use an appropriate PEC or PMC boundary condition in r? In any case this wouldn't work for the Bessel sources.)

However, an even more efficient approach would be to use the principle of equivalence: have a current source surface that completely surrounds the object. Then you should be able to make the computational cell quite small (the only limitation being the inexactness of our PML in the r direction).

[stevengj]

(Eventually, we would like a plot of the error vs. m, and/or the magnitude of the contributions vs m. I would naively expect that it will converge with a power law, so plot on a log–log scale, but it will get worse for increasingly oblique angles.)

[mochen4]

I managed to create the equivalent planewave sources in cylindrical coordinate. Here is the plot of Ep field at m=1, for example:

In comparison, with one long line source of $-\frac1{2}E_p + \frac{i}{2}E_r$ (corresponding to m=1 term of a $E_y$ source), the plot of Ep looks like:

Notice that the colorbar suggests that the amplitude is 0.5 for the equivalent source, but is around 0.3 for the line source. There are also more artifacts for the E fields compared to the H fields; here is a plot of the Hp field, for example:

On the other hand, if I use a line source of H instead, the E field will then look cleaner than the H field.

[stevengj]

It should be exactly a factor of 2 smaller if you had an infinite line source, with just an electric current equal to the planewave amplitude — i.e. if you start with an infinite line equivalent-current source and then set the magnetic currents to zero.

For an oblique planewave, you would:

1. Write the E and H fields of the planewave exp(ikx)
2. Expand each exp(ikx) term in a Jacobi–Anger expansion. Now you have the E and H fields in terms of Bessel functions
3. For each m, use the Bessel E and H fields to construct equivalent sources

(Unlike Cartesian coordinates, I'm not sure that there is any way to get planar wavefronts with just a line source, even if it extends into the PML, except in the limit as the cell radius goes to infinity. The reason is that a finite cylinder with any of the boundary conditions currently in Meep does not support a planewave solution, even along the z direction. One might be able to construct boundary conditions for which this works, however.)

[mochen4]

I created the equivalent sources and the current agreement looks good:

Angle of incidence	$$|m|=0$$	$$|m| \le 1$$	$$|m| \le 2$$	$$|m| \le 3$$	$$|m| \le 4$$	$$|m| \le 5$$
$$0^\circ$$	0.	34.07353489	34.07353489	34.07353489	34.07353489	34.07353489
$$5^\circ$$	1.17112887	33.20208782	33.95312684	33.9596888	33.95971151	33.95971153
$$7.5^\circ$$	2.498413	32.20815769	33.80061922	33.8324966	33.83274769	33.83274833
$$10^\circ$$	4.12241563	30.97946575	33.58485588	33.67994189	33.68129492	33.68130108
$$15^\circ$$	7.49597576	28.34131234	32.96926081	33.37707925	33.39075734	33.39090213

But the errors still increase slowly with angles, so there still might be a small bug somewhere. Out of an abundance of caution, is there any other test I can check?

For reference, here is a plot of the field at $15^\circ$ incidence in air:

[codecov-commenter]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 73.81%. Comparing base (ae9cb64) to head (77e919c).
Report is 46 commits behind head on master.

❗ Your organization needs to install the Codecov GitHub app to enable full functionality.

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #2663      +/-   ##
==========================================
- Coverage   73.81%   73.81%   -0.01%     
==========================================
  Files          18       18              
  Lines        5294     5423     +129     
==========================================
+ Hits         3908     4003      +95     
- Misses       1386     1420      +34     

see 9 files with indirect coverage changes

[stevengj]

Might be worth doubling the resolution to see if that reduces the angle dependence of the results.

[mochen4]

Here are the results at resolution=100 and the agreement looks pretty good with relative error less than 1%

Angle of incidence	$$|m|=0$$	$$|m| \le 1$$	$$|m| \le 2$$	$$|m| \le 3$$	$$|m| \le 4$$	$$|m| \le 5$$
$$0^\circ$$	0.	34.29346804	34.29346804	34.29346804	34.29346804	34.29346804
$$5^\circ$$	1.26307288	33.45160194	34.24236427	34.24941477	34.24943828	34.2494383
$$7.5^\circ$$	2.69179486	32.49031421	34.16626883	34.20051478	34.20077483	34.20077546
$$10^\circ$$	4.43472468	31.29992238	34.03998618	34.1421106	34.1435132	34.14351927
$$15^\circ$$	8.02426432	28.73136299	33.58593652	34.02343398	34.03764361	34.03778666
$$30^\circ$$	10.4188277	24.67896853	30.80532103	33.60631157	34.085187	34.10838805
$$45^\circ$$	8.48250479	21.42266376	28.44485281	32.1072768	33.9610158	34.23488559

[oskooi]

LGTM.