WIP: near2far for cylindrical coordinates. #1090
Conversation
|
To test whether this PR is working correctly or not, the following example computes the far fields at the focal length (known analytically) of a binary-phase zone plate as described in http://zoneplate.lbl.gov/theory. (The specifications of the zone plate are similar to the binary phase grating in Tutorial/Mode Decomposition/Diffraction Spectrum of a Binary Grating.) However, the far fields do not show any focusing behavior. This is unexpected and could either be due to a bug in: (1) the simulation script or (2) For reference, the width of the focal spot is 1.3 um for this configuration computed using eq. 14 of http://zoneplate.lbl.gov/theory. import meep as mp
import numpy as np
import math
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt
resolution = 25 # pixels/um
dpml = 1.0 # PML thickness
dsub = 2.0 # substrate thickness
dpad = 2.0 # padding betweeen zone plate and PML
zh = 0.5 # zone-plate height
zN = 15 # number of zones (odd numbers: \pi phase shift, even numbers: none)
focal_length = 200 # focal length of zone plate
spot_length = 100 # far-field line length
ff_npts = 50 # number of far-field points to compute along spot_length
pml_layers = [mp.PML(thickness=dpml)]
wvl_cen = 0.5 frq_cen = 1/wvl_cen
dfrq = 0.2*frq_cen
## ref: eq. 7 of http://zoneplate.lbl.gov/theory
r = [math.sqrt(n*wvl_cen*(focal_length+n*wvl_cen/4)) for n in range(1,zN+1)]
sr = r[-1]+dpad+dpml
sz = dpml+dsub+zh+dpad+dpml
cell_size = mp.Vector3(sr,0,sz)
sources = [mp.Source(mp.GaussianSource(frq_cen,fwidth=dfrq,is_integrated=True),
component=mp.Er,
center=mp.Vector3(0.5*sr,0,-0.5*sz+dpml),
size=mp.Vector3(sr)),
mp.Source(mp.GaussianSource(frq_cen,fwidth=dfrq,is_integrated=True),
component=mp.Ep,
center=mp.Vector3(0.5*sr,0,-0.5*sz+dpml),
size=mp.Vector3(sr),
amplitude=-1j)]
glass = mp.Medium(index=1.5)
geometry = [mp.Block(material=glass,
size=mp.Vector3(sr,0,dpml+dsub),
center=mp.Vector3(0.5*sr,0,-0.5*sz+0.5*(dpml+dsub)))]
for n in range(zN-1,-1,-1):
geometry.append(mp.Block(material=glass if n % 2 == 0 else mp.vacuum,
size=mp.Vector3(r[n],0,zh),
center=mp.Vector3(0.5*r[n],0,-0.5*sz+dpml+dsub+0.5*zh)))
sim = mp.Simulation(cell_size=cell_size,
boundary_layers=pml_layers,
resolution=resolution,
sources=sources,
geometry=geometry,
dimensions=mp.CYLINDRICAL,
m=-1)
n2f_obj = sim.add_near2far(frq_cen, 0, 1, mp.Near2FarRegion(center=mp.Vector3(0.5*sr,0,0.5*sz-dpml),size=mp.Vector3(sr)))
sim.run(mp.at_beginning(mp.output_epsilon),until_after_sources=50)
ff = sim.get_farfields(n2f_obj, ff_npts/spot_length, center=mp.Vector3(z=-0.5*sz+dpml+dsub+zh+focal_length),size=mp.Vector3(z=spot_length))
energy_density = np.absolute(ff['Ex'])**2+np.absolute(ff['Ey'])**2+np.absolute(ff['Ez'])**2
z = np.linspace(focal_length-0.5*spot_length,focal_length+0.5*spot_length,ff_npts)
plt.figure(dpi=200)
plt.semilogy(z,energy_density,'bo-')
plt.xlabel('z coordinate (um)')
plt.ylabel(r'energy density of far-field electric fields, |E|^2')
plt.tight_layout()
plt.savefig('zone_plate_farfields.png')structure
energy density of far-fields at focal length
|
|
Are you using the latest version of this PR? |
|
Yes: using latest version of this PR with commit message "use meep::pi". As a further diagnostic, the following is the output of the energy density of the electric field separated by component: print(np.absolute(ff['Ex'])**2)
print(np.absolute(ff['Ey'])**2)
print(np.absolute(ff['Ez'])**2)The output contains non-zero values for Ex and Ey which are the same at every point along a line in the focus direction. This is wrong. |
|
I think I've tracked down what the problem is: which produces the following for the above test: This explains why the values of the energy density shown in the previous comment are all identical. Since there does not seem to be anything wrong in the call to |
|
There was a bug in |
|
Working now: the incident planewave is being focused at 200 um as expected.
|
|
Great! Would be good to do a more quantitative test by comparing the near2far predictions with an explicit calculation using a larger cell, the same way that we tested the Cartesian near2far feature. |
|
It would also be nice to know the value of |
|
convergence stats for tol = 1e-3 tol = 1e-5 tol = 1e-12 tol = 1e-16 It's only when the tolerance is reduced to 1e-16 that non-zero values for Nevertheless, in all cases, the far-field profile is identical to the following.
|
Probably we should change the starting What do the "(99.95%)" numbers in your output indicate? What is 57060? |
|
The percentage values refer to the fraction of all converged calls of As an alternative to doubling |
Note that the re-using of points from the previous It's probably not worth it. Let's just change |
|
As additional validation that this PR is working correctly, the far fields obtained using an actual enlarged cell are proportional (by roughly a constant factor) to those obtained using
import meep as mp
import numpy as np
import math
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt
resolution = 25 # pixels/um
dpml = 1.0 # PML thickness
dsub = 2.0 # substrate thickness
dpad = 2.0 # padding betweeen zone plate and PML
zh = 0.5 # zone-plate height
zN = 15 # number of zones (odd numbers: \pi phase shift, even numbers: none)
focal_length = 200 # focal length of zone plate
spot_length = 100 # far-field line length
pml_layers = [mp.PML(thickness=dpml)]
wvl_cen = 0.5
frq_cen = 1/wvl_cen
dfrq = 0.2*frq_cen
## ref: eq. 7 of http://zoneplate.lbl.gov/theory
r = [math.sqrt(n*wvl_cen*(focal_length+n*wvl_cen/4)) for n in range(1,zN+1)]
sr = r[-1]+dpad+dpml
sz = dpml+dsub+zh+focal_length+0.5*spot_length+dpad+dpml
cell_size = mp.Vector3(sr,0,sz)
sources = [mp.Source(mp.GaussianSource(frq_cen,fwidth=dfrq,is_integrated=True),
component=mp.Er,
center=mp.Vector3(0.5*sr,0,-0.5*sz+dpml),
size=mp.Vector3(sr)),
mp.Source(mp.GaussianSource(frq_cen,fwidth=dfrq,is_integrated=True),
component=mp.Ep,
center=mp.Vector3(0.5*sr,0,-0.5*sz+dpml),
size=mp.Vector3(sr),
amplitude=-1j)]
glass = mp.Medium(index=1.5)
geometry = [mp.Block(material=glass,
size=mp.Vector3(sr,0,dpml+dsub),
center=mp.Vector3(0.5*sr,0,-0.5*sz+0.5*(dpml+dsub)))]
for n in range(zN-1,-1,-1):
geometry.append(mp.Block(material=glass if n % 2 == 0 else mp.vacuum,
size=mp.Vector3(r[n],0,zh),
center=mp.Vector3(0.5*r[n],0,-0.5*sz+dpml+dsub+0.5*zh)))
sim = mp.Simulation(cell_size=cell_size,
boundary_layers=pml_layers,
resolution=resolution,
sources=sources,
geometry=geometry,
dimensions=mp.CYLINDRICAL,
m=-1)
dft_obj = sim.add_dft_fields([mp.Er,mp.Ep,mp.Ez], frq_cen, frq_cen, 1, center=mp.Vector3(z=-0.5*sz+dpml+dsub+zh+focal_length), size=mp.Vector3(z=spot_length))
sim.run(until_after_sources=400)
energy_density = np.absolute(sim.get_dft_array(dft_obj,mp.Er,0))**2+np.absolute(sim.get_dft_array(dft_obj,mp.Ep,0))**2+np.absolute(sim.get_dft_array(dft_obj,mp.Ez,0))**2
z = np.linspace(focal_length-0.5*spot_length,focal_length+0.5*spot_length,len(energy_density))
if mp.am_master():
plt.figure(dpi=200)
plt.semilogy(z,energy_density,'bo-')
plt.xlabel('z coordinate (um)')
plt.ylabel(r'energy density of far-field electric fields, |E|^2')
plt.title("focusing peropties of a binary-phase zone plate (extended cell)")
plt.tight_layout()
plt.savefig("zone_plate_extended_cell_farfields_zn{}.png".format(zN)) |
|
If they are only proportional, not nearly identical, then something is wrong. What is the proportionality factor? |
|
The far fields from the |
No, that shouldn't cause any difference bigger than discretization errors. |
|
I noticed that a factor of 2π was being included twice, so in the energy you would have gotten an extra factor of (2π)². Can you try it again with the updated PR? |
|
There is a bug in commit |
|
For testing purposes, I would use a much smaller cell, e.g. predicting the fields 10µm away. Not only will this be a lot faster, but it should also be a lot easier to get good agreement — when you propagate for a really long distance in a brute-force calculations, numerical dispersion effects become larger and larger, so you will probably need higher resolution to get agreement for a large cell. |
|
Even with a smaller cell (focal length of 10 µm rather than 200 µm), the
|
|
Can you try a point with r > 0? |
|
The discrepancy between the two methods becomes larger for r > 0 as shown below for r=0.2 and r=4.6.
|
|
Okay, so it's not a constant factor at all, it's some more nontrivial bug. 😿 |
|
Not sure what the problem here is. One thing to try would be an You could also try to increase (The basic problem here is that there are a lot of sign errors and things that would still give qualitatively reasonable results.) |
|
Note: as the point where you are trying to compute the fields gets closer to the near2far surface, the computation time will get larger because the required There are many ways to mitigate this, and in fact there is a whole literature on this sort of integral because it arises for integral-equation methods (like Homer's SCUFF code), but they will all be much more complicated than what we are doing now. I doubt that it is worth it because — once we get this debugged — the practical application of the near2far code is to calculate the field at far-away points, not at nearby points. |
|
The test involving The same is true for the magnetic fields: Hr is The brute-force DFT fields are as expected: summary: in the near2far calculation, the far fields are not being comptued correctly for the case of addendum: Eφ and Hr at
import meep as mp
import numpy as np
import math
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt
## run 1: near2far calculation
resolution = 25 # pixels/um
dpml = 1.0 # PML thickness
pml_layers = [mp.PML(thickness=dpml)]
wvl_cen = 0.5
frq_cen = 1/wvl_cen
r = 5
z = 8
sr = r+dpml
sz = dpml+z+dpml
cell_size = mp.Vector3(sr,0,sz)
sources = [mp.Source(mp.GaussianSource(frq_cen,fwidth=0.2*frq_cen),
component=mp.Ep,
center=mp.Vector3(0.5*(sr-dpml),0,-0.5*sz+dpml),
size=mp.Vector3())]
sim = mp.Simulation(cell_size=cell_size,
boundary_layers=pml_layers,
resolution=resolution,
sources=sources,
dimensions=mp.CYLINDRICAL,
m=0)
n2f_obj = sim.add_near2far(frq_cen, 0, 1, mp.Near2FarRegion(center=mp.Vector3(0.5*(sr-dpml),0,0.5*sz-dpml),size=mp.Vector3(sr-dpml)))
sim.run(until_after_sources=100)
ff_distance = 20
ff = sim.get_farfields(n2f_obj, resolution, center=mp.Vector3(sr-dpml-0.5*(sr-dpml),0,-0.5*sz+dpml+z+ff_distance), size=mp.Vector3(sr-dpml))
E2_n2f = np.absolute(ff['Ex'])**2+np.absolute(ff['Ey'])**2+np.absolute(ff['Ez'])**2
sim.reset_meep()
## run 2: brute-force calculation
sz = dpml+z+ff_distance+dpml
cell_size = mp.Vector3(sr,0,sz)
sources = [mp.Source(mp.GaussianSource(frq_cen,fwidth=0.2*frq_cen),
component=mp.Ep,
center=mp.Vector3(0.5*(sr-dpml),0,-0.5*sz+dpml),
size=mp.Vector3())]
sim = mp.Simulation(cell_size=cell_size,
boundary_layers=pml_layers,
resolution=resolution,
sources=sources,
dimensions=mp.CYLINDRICAL,
m=0)
dft_obj = sim.add_dft_fields([mp.Er,mp.Ep,mp.Ez], frq_cen, frq_cen, 1, center=mp.Vector3(0.5*(sr-dpml),0,-0.5*sz+dpml+z+ff_distance), size=mp.Vector3(sr-dpml))
sim.run(until_after_sources=200)
E2_dft = np.absolute(sim.get_dft_array(dft_obj,mp.Er,0))**2+np.absolute(sim.get_dft_array(dft_obj,mp.Ep,0))**2+np.absolute(sim.get_dft_array(dft_obj,mp.Ez,0))**2
if mp.am_master():
plt.figure(dpi=200)
plt.semilogy(np.linspace(0,sr-dpml,len(E2_n2f)),E2_n2f,'bo-',label="near2far")
plt.semilogy(np.linspace(0,sr-dpml,len(E2_dft)),E2_dft,'ro-',label="extended cell")
plt.grid(True,axis="y",which="both",ls="-")
plt.legend(loc='upper right')
plt.xlabel('r coordinate (μm)')
plt.ylabel(r'energy density of far-field electric fields, |E|$^2$')
plt.title("empty cell in cylindrical coordinates with m=0")
plt.tight_layout()
plt.savefig("empty_cyl_compare_farfields_res{}_fd{}.png".format(resolution,ff_distance)) |
|
The NaN is arising just from The issue is that the DFT weight includes the |
|
The next thing we might try is to test the
Ideally they should be the same (up to an overall normalization). |
|
The following are the results for the Green's function for ε=2 everywhere computed in C++ using (1) the brute-force approach via the CW solver and (2) directly from One thing to note:
#include <stdio.h>
#include <stdlib.h>
#include <meep.hpp>
using namespace meep;
using std::complex;
using std::polar;
double two(const vec &) { return 2.0; }
int check_greencyl(double sr, double sz, double res, component c0, double m, double tol) {
const double dpml = 2.0;
grid_volume gv = volcyl(sr+dpml, dpml+sz+dpml, res);
gv.center_origin();
structure s(gv, two, pml(dpml));
fields f(&s, m);
double w = 0.5;
continuous_src_time src(w);
vec x0 = veccyl(0.5*sr,-0.5*sz);
f.add_point_source(c0, src, x0);
f.solve_cw(1e-6);
f.output_hdf5(c0, gv.surroundings());
const int N = 10;
double dr = sr / N;
complex<double> F[N], F0[N], EH[6];
complex<double> phase = polar(1.0, (4 * w * f.dt) * pi);
for (int i = 0; i < N; ++i) {
double rr = dr * i;
vec x = veccyl(rr, 0.5*sz-0.1);
F[i] = f.get_field(Ep, x) * phase;
greencyl(EH, x, w, 2.0, 1.0, x0, c0, 1.0, m, tol);
F0[i] = EH[1]; // Ey = Ep for \phi = 0 (rz = xz) plane
master_printf("green:, %0.2f, %0.10f, %0.10f, %0.10f\n",rr,abs(F[i]),abs(F0[i]),abs(F[i])/abs(F0[i]));
}
return 0;
}
int main(int argc, char **argv) {
initialize mpi(argc, argv);
check_greencyl(5.0, 10.0, 40.0, Ep, 0, 1e-6);
return 0;
}output
|
|
Four additional notes for debugging:
|
|
I can confirm that the fields it's returning for an #include "meep_internals.hpp"
#include <complex>
#include <stdio.h>
using namespace meep;
using namespace std;
int main(void) {
std::complex<double> EH[6];
vec x0 = veccyl(1,2), x = veccyl(2,5);
double freq = 0.5;
greencyl(EH, x, freq, 1.0, 1.0, x0, Ep, 1.0, 0.0, 1e-6);
printf("Er=%g%+gi, Ep=%g%+gi, Ez=%g%+gi\n",
real(EH[0]), imag(EH[0]),
real(EH[1]), imag(EH[1]),
real(EH[2]), imag(EH[2]));
return 0;
} |
…t them some day, and to make testing easier
|
Just pushed a fix to the symmetry — was a simple bug where I was evaluating at the wrong angles. |
|
Seems to be working now: the ratio of the magnitude of |
|
Great! I'm expecting a constant factor in this test because the normalizations are different. The |
|
Here is the output of the real and imaginary parts of the ratio obtained using: master_printf("green:, %0.2f, %g%+gi, %g%+gi, %g%+gi\n",rr,real(F[i]),imag(F[i]),real(F0[i]),imag(F0[i]),real(F[i]/F0[i]),imag(F[i]/F0[i])); |
|
In particular, the (Your imaginary parts look small, around 2% of the real parts, so I think that is just a discretization error in the phase, which is expected.) |
|
Changing to |
|
As an additional check, the m=1 case is also working now (even though the CW convergence failed): |
|
The far fields at the focal length of the binary-phase zone plate example are also now showing good agreement for a monitor line along both r and z as shown below.
|
|
I can merge this now as an undocumented new feature if you want, and then you can make a separate PR with a test and documentation? |
|
Yes, we can merge this now. |
* initial stab at a greencyl function * whoops * formatting * use meep::pi * fix get_farfields_array in cylindrical coords * slight code cleanup * rm redundant 2pi factor * smaller N0 * add some comments to greencyl * typo in comment * whoops, 2/N is integer division in C * add green functions to meep internals header in case we want to export them some day, and to make testing easier * whoops, fix evaluation angles * clarify normalization of greencyl
This is a rough draft to fix #1086 using the simplest possible method: integrating
green3dnumerically over φ to get the corresponding cylindrical function.(Integrates by repeatedly doubling the number of φ points to a tolerance of
1e-3or 2^16 points, whichever comes first, with a simple Euler/trapezoidal rule — which should be exponentially convergent for smooth cylindrical integrals like this.)Untested except to check that it compiles, intended as a first draft.