WIP: periodic green's functions #769
Merged
Conversation
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
|
There is likely a simple bug in this PR because it does not appear to actually do anything. This is demonstrated using the script below which computes the far-field intensity profile relative to vacuum at λ=0.5 μm along a line for a binary grating. The test involves verifying that the diffraction spectra contains only odd orders, m=±1 (±2.9°), ±3 (±8.6°), ±5 (±14.5°), etc., and that the ratio of the peaks for any two orders m1 and m2 matches the analytic result (m2)2/(m1)2. import meep as mp
import math
import matplotlib.pyplot as plt
import numpy as np
resolution = 20 # pixels/μm
dpml = 1.0 # PML thickness
dsub = 3.0 # substrate thickness
dpad = 3.0 # padding between grating and PML
gp = 10.0 # grating period
gh = 0.5 # grating height
gdc = 0.5 # grating duty cycle
sx = dpml+dsub+gh+dpad+dpml
pml_layers = [mp.PML(thickness=dpml,direction=mp.X)]
wvl = 0.5 # center wavelength
fcen = 1/wvl # center frequency
df = 0.05*fcen # frequency width
src_pt = mp.Vector3(-0.5*sx+dpml+0.5*dsub)
sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df), component=mp.Ez, center=src_pt)]
k_point = mp.Vector3(0,0,0)
glass = mp.Medium(index=1.5)
sim = mp.Simulation(resolution=resolution,
cell_size=mp.Vector3(sx),
boundary_layers=pml_layers,
k_point=k_point,
default_material=glass,
sources=sources)
mon_pt = mp.Vector3(0.5*sx-dpml-0.5*dpad)
n2f_mon = sim.add_near2far(fcen, 0, 1, mp.Near2FarRegion(center=mon_pt, size=mp.Vector3()))
sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, mon_pt, 1e-9))
ff_x = 1e8 # location of line (10 m away from grating)
ff_angle = 20 # half angle subtended by the line
ff_y = 2*ff_x*math.tan(math.radians(ff_angle)) # length of line
ff_npts = 800 # number of evenly-spaced points along line to compute far fields
ff_res = ff_npts/ff_y
ff_source = sim.get_farfields(n2f_mon, ff_res, center=mp.Vector3(ff_x), size=mp.Vector3(y=ff_y))
sim.reset_meep()
num_cells = 15 # number of unit cells in super cell
sy = num_cells*gp
sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df), component=mp.Ez, center=src_pt, size=mp.Vector3(y=sy))]
geometry = [mp.Block(material=glass, size=mp.Vector3(dpml+dsub,mp.inf,mp.inf), center=mp.Vector3(-0.5*sx+0.5*(dpml+dsub)))]
for n in range(num_cells):
geometry.append(mp.Block(material=glass, size=mp.Vector3(gh,gdc*gp,mp.inf), center=mp.Vector3(-0.5*sx+dpml+dsub+0.5*gh,-0.5*sy+(n+0.5)*gp)))
symmetries=[mp.Mirror(mp.Y)]
sim = mp.Simulation(resolution=resolution,
cell_size=mp.Vector3(sx,sy),
boundary_layers=pml_layers,
geometry=geometry,
k_point=k_point,
sources=sources,
symmetries=symmetries)
n2f_mon = sim.add_near2far(fcen, 0, 1, mp.Near2FarRegion(center=mon_pt, size=mp.Vector3(y=sy)))
sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, mon_pt, 1e-9))
ff_grating = sim.get_farfields(n2f_mon, ff_res, center=mp.Vector3(ff_x), size=mp.Vector3(y=ff_y))
ff_ys = np.linspace(-0.5*ff_y,0.5*ff_y,ff_npts)
ang = np.degrees(np.arctan(ff_ys/ff_x))
ff_enh = abs(ff_grating['Ez'][0])**2/abs(ff_source['Ez'][0])**2
plt.figure()
plt.plot(ang,ff_enh,'bo-')
plt.xlabel('angle (degrees)')
plt.ylabel('far-field enhancement (relative to vacuum)');
plt.title('number of unit cells: {}'.format(num_cells))
plt.show()In the first part of the test, which does not involve this PR, the unit cell is periodically tiled to create a supercell. The results are shown below and demonstrate that as the number of unit cells is increased (
In the second test which is for this PR, only a single unit cell with periodic boundaries is used. The far-field spectra is shown below. It is identical to the left plot from the figure above. This indicates that this PR is having no effect.
The correct result would be for the far-field spectra computed for a single unit cell using this PR to exactly match that of a 20 unit-cell supercell without this PR. |
oskooi
reviewed
Mar 16, 2019
src/near2far.cpp
Outdated
| for (int i = 0; i < 2; ++i) { | ||
| if (has_direction(v.dim, fd[i]) && | ||
| boundaries[High][fd[i]] == Periodic && boundaries[Low][fd[i]] == Periodic && | ||
| float(w->v.in_direction(fd[i])) == float(v.in_direction(fd[i]))) { |
There was a problem hiding this comment.
It seems the problem lies in this particular if statement, particularly the third and final test condition involving in_direction, which always evaluates to false and thus the subsequent initialization of the arrays periodic_d, periodic_n, period, and periodic_k to non-zero values is never performed.
stevengj
commented
Mar 21, 2019
src/near2far.cpp
Outdated
| xs.set_direction(periodic_d[1], x0.in_direction(periodic_d[1]) + i1*period[1]); | ||
| double phase = phase0 + i1*periodic_k[1]; | ||
| std::complex<double> cphase = std::polar(1.0, phase); | ||
| green(EH6, x, freq, eps, mu, x0, c0, f->dft[Nfreq * idx_dft + i]); |
There was a problem hiding this comment.
whoops: this should be xs, not x0!
|
With the one fix from
|
|
For reference, the following script was used to generate results for the unit-cell calculation. Note that import meep as mp
import math
import matplotlib.pyplot as plt
import numpy as np
resolution = 20 # pixels/μm
dpml = 1.0 # PML thickness
dsub = 3.0 # substrate thickness
dpad = 3.0 # padding between grating and PML
gp = 10.0 # grating period
gh = 0.5 # grating height
gdc = 0.5 # grating duty cycle
sx = dpml+dsub+gh+dpad+dpml
sy = gp
pml_layers = [mp.PML(thickness=dpml,direction=mp.X)]
wvl = 0.5 # center wavelength
fcen = 1/wvl # center frequency
df = 0.05*fcen # frequency width
src_pt = mp.Vector3(-0.5*sx+dpml+0.5*dsub)
sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df), component=mp.Ez, center=src_pt)]
k_point = mp.Vector3(0,0,0)
glass = mp.Medium(index=1.5)
sim = mp.Simulation(resolution=resolution,
cell_size=mp.Vector3(sx,0,0),
boundary_layers=pml_layers,
k_point=k_point,
default_material=glass,
sources=sources)
mon_pt = mp.Vector3(0.5*sx-dpml-0.5*dpad)
n2f_mon = sim.add_near2far(fcen, 0, 1, mp.Near2FarRegion(center=mon_pt, size=mp.Vector3()))
sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, mon_pt, 1e-9))
ff_x = 1e8
ff_angle = 20
ff_y = 2*ff_x*math.tan(math.radians(ff_angle)) # first 9 modes
ff_npts = 800
ff_res = ff_npts/ff_y
ff_source = sim.get_farfields(n2f_mon, ff_res, center=mp.Vector3(ff_x), size=mp.Vector3(y=ff_y))
sim.reset_meep()
sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df), component=mp.Ez, center=src_pt, size=mp.Vector3(y=sy))]
geometry = [mp.Block(material=glass, size=mp.Vector3(dpml+dsub,mp.inf,mp.inf), center=mp.Vector3(-0.5*sx+0.5*(dpml+dsub))),
mp.Block(material=glass, size=mp.Vector3(gh,gdc*gp,mp.inf), center=mp.Vector3(-0.5*sx+dpml+dsub+0.5*gh))]
sim = mp.Simulation(resolution=resolution,
cell_size=mp.Vector3(sx,sy,0),
boundary_layers=pml_layers,
geometry=geometry,
k_point=k_point,
sources=sources)
n2f_mon = sim.add_near2far(fcen, 0, 1, mp.Near2FarRegion(center=mon_pt, size=mp.Vector3(y=sy)))
sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, mon_pt, 1e-9))
ff_grating = sim.get_farfields(n2f_mon, ff_res, center=mp.Vector3(ff_x), size=mp.Vector3(y=ff_y))
ff_ys = np.linspace(-0.5*ff_y,0.5*ff_y,ff_npts)
ang = np.degrees(np.arctan(ff_ys/ff_x))
plt.figure()
plt.plot(ang,abs(ff_grating['Ez'][0])**2/abs(ff_source['Ez'][0])**2,'bo-')
plt.xlabel('angle (degrees)')
plt.ylabel('far-field enhancement (relative to source)');
plt.title('unit cell with periodic boundaries using PR #769')
plt.show() |
|
To use it, pass |
|
|
|
updated script for unit-cell calculation which includes import meep as mp
import math
import matplotlib.pyplot as plt
import numpy as np
resolution = 20 # pixels/μm
dpml = 1.0 # PML thickness
dsub = 3.0 # substrate thickness
dpad = 3.0 # padding between grating and PML
gp = 10.0 # grating period
gh = 0.5 # grating height
gdc = 0.5 # grating duty cycle
sx = dpml+dsub+gh+dpad+dpml
sy = gp
pml_layers = [mp.PML(thickness=dpml,direction=mp.X)]
wvl = 0.5 # center wavelength
fcen = 1/wvl # center frequency
df = 0.05*fcen # frequency width
src_pt = mp.Vector3(-0.5*sx+dpml+0.5*dsub)
sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df), component=mp.Ez, center=src_pt)]
k_point = mp.Vector3(0,0,0)
glass = mp.Medium(index=1.5)
sim = mp.Simulation(resolution=resolution,
cell_size=mp.Vector3(sx),
boundary_layers=pml_layers,
k_point=k_point,
default_material=glass,
sources=sources)
mon_pt = mp.Vector3(0.5*sx-dpml-0.5*dpad)
n2f_mon = sim.add_near2far(fcen, 0, 1, mp.Near2FarRegion(center=mon_pt, size=mp.Vector3()))
sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, mon_pt, 1e-9))
ff_x = 1e8
ff_angle = 21
ff_y = 2*ff_x*math.tan(math.radians(ff_angle)) # first 9 modes
ff_npts = 800
ff_res = ff_npts/ff_y
ff_source = sim.get_farfields(n2f_mon, ff_res, center=mp.Vector3(ff_x), size=mp.Vector3(y=ff_y))
sim.reset_meep()
sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df), component=mp.Ez, center=src_pt, size=mp.Vector3(y=sy))]
geometry = [mp.Block(material=glass, size=mp.Vector3(dpml+dsub,mp.inf,mp.inf), center=mp.Vector3(-0.5*sx+0.5*(dpml+dsub))),
mp.Block(material=glass, size=mp.Vector3(gh,gdc*gp,mp.inf), center=mp.Vector3(-0.5*sx+dpml+dsub+0.5*gh))]
# symmetries=[mp.Mirror(mp.Y)]
symmetries = []
sim = mp.Simulation(resolution=resolution,
cell_size=mp.Vector3(sx,sy),
boundary_layers=pml_layers,
geometry=geometry,
k_point=k_point,
sources=sources,
symmetries=symmetries)
nd = 10
n2f_mon = sim.add_near2far(fcen, 0, 1, mp.Near2FarRegion(center=mon_pt, size=mp.Vector3(y=sy)), nperiods=nd)
sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, mon_pt, 1e-9))
ff_grating = sim.get_farfields(n2f_mon, ff_res, center=mp.Vector3(ff_x), size=mp.Vector3(y=ff_y))
ff_ys = np.linspace(-0.5*ff_y,0.5*ff_y,ff_npts)
ang = np.degrees(np.arctan(ff_ys/ff_x))
plt.figure()
plt.plot(ang,abs(ff_grating['Ez'][0])**2/abs(ff_source['Ez'][0])**2,'bo-')
plt.xlabel('angle (degrees)')
plt.ylabel('far-field enhancement (relative to source)');
plt.title('unit cell with periodic boundaries, nperiods = 7, using PR #769 \n without symmetries')
plt.show() |
|
Working now: results are identical with and without
|
This was referenced Mar 26, 2019
bencbartlett
pushed a commit
to bencbartlett/meep
that referenced
this pull request
Sep 9, 2021
* first attempt to hack in periodic green's functions * typo * make nperiods an option for near2far * only do lattice sum in non-empty dirs * simplify kw handling in add_near2far * fix symmetry again
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
First attempt (currently untested) to hack in support for periodic Green's functions (closes #763).
That is, it detects whether any directions of the near-to-far surface coincide with a periodic direction of the computational cell. If so, in that direction it just does a naive Bloch-periodic summation of the Green's function. Currently, it just sums 10 copies. For far-away surfaces we may need many more, and in principle should do something like Ewald summation, but for now this is hopefully good enough to experiment with.
Update: modulo the bugfix reported below, @oskooi reports that it seems to be working in a quick test. To do:
x0typo noted below.