Closed mochen4 closed 4 months ago
Ensuring stability(#2644) and forcing complex fields (inspired by https://github.com/NanoComp/meep/pull/2726#issuecomment-1899803540) 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:
Good to see this is finally working as expected. A few comments and questions:
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?m
of M = 4
for each angle? should M
be a constant independent of the angle of the planewave?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.)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
?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
ofM = 4
for each angle? shouldM
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
toM
? why not sum from0
toM
where the results for the flux for them > 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 theCourant
factor need to be a constant? What about making its value a function ofm
?
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.
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).
(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.)
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.
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:
(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.)
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:
All modified and coverable lines are covered by tests :white_check_mark:
Project coverage is 73.81%. Comparing base (
ae9cb64
) to head (77e919c
). Report is 46 commits behind head on master.
:exclamation: Your organization needs to install the Codecov GitHub app to enable full functionality.
Might be worth doubling the resolution to see if that reduces the angle dependence of the results.
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 |
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![oblique_cyl](https://github.com/NanoComp/meep/assets/25192039/a948a005-a9e8-4656-8b02-c58715d533d3)
Part of the source around $r=0$ is currently excluded for $|m|>1$ due to #2644. In addition, as noted in https://github.com/NanoComp/meep/pull/2538#issuecomment-1686430983, 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|
atresolution=100
.A higher cutoff should have better convergence, but as mentioned in https://github.com/NanoComp/meep/issues/2644#issuecomment-1743208682, a higher cutoff currently required excluding more pixels around $r=0$.