From 15b3c567cbd297aa4358994d28fa67db0a117573 Mon Sep 17 00:00:00 2001 From: Ardavan Oskooi Date: Wed, 16 Mar 2022 19:02:55 -0700 Subject: [PATCH] add support for out-of-plane wavevector in 2d cell to fields::get_eigenmode (#1968) * add support for out-of-plane wavevector in 2d cell to fields::get_eigenmode * add unit test and apply phase fix only for real fields and beta * sum S- and P-polarized diffraction orders when kz is non-zero --- python/tests/test_binary_grating.py | 256 ++++++++++++++++++++++------ python/tests/test_special_kz.py | 41 +++-- src/mpb.cpp | 33 +++- 3 files changed, 263 insertions(+), 67 deletions(-) diff --git a/python/tests/test_binary_grating.py b/python/tests/test_binary_grating.py index 380c719e9..51b5a0e69 100644 --- a/python/tests/test_binary_grating.py +++ b/python/tests/test_binary_grating.py @@ -1,47 +1,61 @@ +# -*- coding: utf-8 -*- + import unittest +import parameterized import meep as mp import math import cmath import numpy as np + class TestEigCoeffs(unittest.TestCase): + @classmethod + def setUpClass(cls): + cls.resolution = 30 # pixels/μm - def run_binary_grating_oblique(self, theta): + cls.dpml = 1.0 # PML thickness + cls.dsub = 1.0 # substrate thickness + cls.dpad = 1.0 # padding thickness between grating and PML + cls.gp = 6.0 # grating period + cls.gh = 0.5 # grating height + cls.gdc = 0.5 # grating duty cycle - resolution = 30 # pixels/um + cls.sx = cls.dpml+cls.dsub+cls.gh+cls.dpad+cls.dpml + cls.sy = cls.gp - dpml = 1.0 # PML thickness - dsub = 1.0 # substrate thickness - dpad = 1.0 # length of padding between grating and pml - gp = 6.0 # grating period - gh = 0.5 # grating height - gdc = 0.5 # grating duty cycle + cls.cell_size = mp.Vector3(cls.sx,cls.sy,0) - sx = dpml+dsub+gh+dpad+dpml - sy = gp + # replace anisotropic PML with isotropic Absorber to + # attenuate parallel-directed fields of oblique source + cls.abs_layers = [mp.Absorber(thickness=cls.dpml,direction=mp.X)] - cell_size = mp.Vector3(sx,sy,0) + wvl = 0.5 # center wavelength + cls.fcen = 1/wvl # center frequency + cls.df = 0.05*cls.fcen # frequency width - # replace anisotropic PML with isotropic Absorber to attenuate parallel-directed fields of oblique source - abs_layers = [mp.Absorber(thickness=dpml,direction=mp.X)] + cls.ng = 1.5 + cls.glass = mp.Medium(index=cls.ng) - wvl = 0.5 # center wavelength - fcen = 1/wvl # center frequency - df = 0.05*fcen # frequency width + cls.geometry = [mp.Block(material=cls.glass, + size=mp.Vector3(cls.dpml+cls.dsub,mp.inf,mp.inf), + center=mp.Vector3(-0.5*cls.sx+0.5*(cls.dpml+cls.dsub),0,0)), + mp.Block(material=cls.glass, + size=mp.Vector3(cls.gh,cls.gdc*cls.gp,mp.inf), + center=mp.Vector3(-0.5*cls.sx+cls.dpml+cls.dsub+0.5*cls.gh,0,0))] - ng = 1.5 - glass = mp.Medium(index=ng) - # rotation angle of incident planewave; CCW about Y axis, 0 degrees along +X axis + @parameterized.parameterized.expand([(0,), (10.7,)]) + def test_binary_grating_oblique(self, theta): + # rotation angle of incident planewave + # counterclockwise (CCW) about Z axis, 0 degrees along +X axis theta_in = math.radians(theta) # k (in source medium) with correct length (plane of incidence: XY) - k = mp.Vector3(math.cos(theta_in),math.sin(theta_in),0).scale(fcen*ng) + k = mp.Vector3(self.fcen*self.ng).rotate(mp.Vector3(0,0,1), theta_in) symmetries = [] eig_parity = mp.ODD_Z if theta_in == 0: - k = mp.Vector3(0,0,0) symmetries = [mp.Mirror(mp.Y)] eig_parity += mp.EVEN_Y @@ -50,68 +64,84 @@ def _pw_amp(x): return cmath.exp(1j*2*math.pi*k.dot(x+x0)) return _pw_amp - src_pt = mp.Vector3(-0.5*sx+dpml+0.3*dsub,0,0) - sources = [mp.Source(mp.GaussianSource(fcen,fwidth=df), - component=mp.Ez, + src_pt = mp.Vector3(-0.5*self.sx+self.dpml,0,0) + sources = [mp.Source(mp.GaussianSource(self.fcen,fwidth=self.df), + component=mp.Ez, # S polarization center=src_pt, - size=mp.Vector3(0,sy,0), + size=mp.Vector3(0,self.sy,0), amp_func=pw_amp(k,src_pt))] - sim = mp.Simulation(resolution=resolution, - cell_size=cell_size, - boundary_layers=abs_layers, + sim = mp.Simulation(resolution=self.resolution, + cell_size=self.cell_size, + boundary_layers=self.abs_layers, k_point=k, - default_material=glass, + default_material=self.glass, sources=sources, symmetries=symmetries) - refl_pt = mp.Vector3(-0.5*sx+dpml+0.5*dsub,0,0) - refl_flux = sim.add_flux(fcen, 0, 1, mp.FluxRegion(center=refl_pt, size=mp.Vector3(0,sy,0))) + refl_pt = mp.Vector3(-0.5*self.sx+self.dpml+0.5*self.dsub,0,0) + refl_flux = sim.add_flux(self.fcen, + 0, + 1, + mp.FluxRegion(center=refl_pt, + size=mp.Vector3(0,self.sy,0))) - sim.run(until_after_sources=100) + sim.run(until_after_sources=mp.stop_when_dft_decayed()) input_flux = mp.get_fluxes(refl_flux) input_flux_data = sim.get_flux_data(refl_flux) sim.reset_meep() - geometry = [mp.Block(material=glass, size=mp.Vector3(dpml+dsub,mp.inf,mp.inf), center=mp.Vector3(-0.5*sx+0.5*(dpml+dsub),0,0)), - mp.Block(material=glass, size=mp.Vector3(gh,gdc*gp,mp.inf), center=mp.Vector3(-0.5*sx+dpml+dsub+0.5*gh,0,0))] - - sim = mp.Simulation(resolution=resolution, - cell_size=cell_size, - boundary_layers=abs_layers, - geometry=geometry, + sim = mp.Simulation(resolution=self.resolution, + cell_size=self.cell_size, + boundary_layers=self.abs_layers, + geometry=self.geometry, k_point=k, sources=sources, symmetries=symmetries) - refl_flux = sim.add_flux(fcen, 0, 1, mp.FluxRegion(center=refl_pt, size=mp.Vector3(0,sy,0))) + refl_flux = sim.add_flux(self.fcen, + 0, + 1, + mp.FluxRegion(center=refl_pt, + size=mp.Vector3(0,self.sy,0))) + sim.load_minus_flux_data(refl_flux,input_flux_data) - tran_pt = mp.Vector3(0.5*sx-dpml-0.5*dpad,0,0) - tran_flux = sim.add_flux(fcen, 0, 1, mp.FluxRegion(center=tran_pt, size=mp.Vector3(0,sy,0))) + tran_pt = mp.Vector3(0.5*self.sx-self.dpml-0.5*self.dpad,0,0) + tran_flux = sim.add_flux(self.fcen, + 0, + 1, + mp.FluxRegion(center=tran_pt, + size=mp.Vector3(0,self.sy,0))) - sim.run(until_after_sources=100) + sim.run(until_after_sources=mp.stop_when_dft_decayed()) - nm_r = np.floor((fcen*ng-k.y)*gp)-np.ceil((-fcen*ng-k.y)*gp) # number of reflected orders + # number of reflected orders + nm_r = np.floor((self.fcen*self.ng-k.y)*self.gp)-np.ceil((-self.fcen*self.ng-k.y)*self.gp) if theta_in == 0: nm_r = nm_r/2 # since eig_parity removes degeneracy in y-direction nm_r = int(nm_r) - res = sim.get_eigenmode_coefficients(refl_flux, range(1,nm_r+1), eig_parity=eig_parity) + res = sim.get_eigenmode_coefficients(refl_flux, + range(1,nm_r+1), + eig_parity=eig_parity) r_coeffs = res.alpha Rsum = 0 for nm in range(nm_r): Rsum += abs(r_coeffs[nm,0,1])**2/input_flux[0] - nm_t = np.floor((fcen-k.y)*gp)-np.ceil((-fcen-k.y)*gp) # number of transmitted orders + # number of transmitted orders + nm_t = np.floor((self.fcen-k.y)*self.gp)-np.ceil((-self.fcen-k.y)*self.gp) if theta_in == 0: nm_t = nm_t/2 # since eig_parity removes degeneracy in y-direction nm_t = int(nm_t) - res = sim.get_eigenmode_coefficients(tran_flux, range(1,nm_t+1), eig_parity=eig_parity) + res = sim.get_eigenmode_coefficients(tran_flux, + range(1,nm_t+1), + eig_parity=eig_parity) t_coeffs = res.alpha Tsum = 0 @@ -123,12 +153,138 @@ def _pw_amp(x): Rflux = -r_flux[0]/input_flux[0] Tflux = t_flux[0]/input_flux[0] + print("refl:, {}, {}".format(Rsum,Rflux)) + print("tran:, {}, {}".format(Tsum,Tflux)) + print("sum:, {}, {}".format(Rsum+Tsum,Rflux+Tflux)) + + self.assertAlmostEqual(Rsum,Rflux,places=2) + self.assertAlmostEqual(Tsum,Tflux,places=2) + self.assertAlmostEqual(Rsum+Tsum,1.00,places=2) + + + @parameterized.parameterized.expand([(13.2,"real/imag"), + (17.7,"complex"), + (21.2, "3d")]) + def test_binary_grating_special_kz(self, theta, kz_2d): + # rotation angle of incident planewave + # counterclockwise (CCW) about Y axis, 0 degrees along +X axis + theta_in = math.radians(theta) + + # k (in source medium) with correct length (plane of incidence: XZ) + k = mp.Vector3(self.fcen*self.ng).rotate(mp.Vector3(0,1,0), theta_in) + + symmetries = [mp.Mirror(mp.Y)] + + def pw_amp(k,x0): + def _pw_amp(x): + return cmath.exp(1j*2*math.pi*k.dot(x+x0)) + return _pw_amp + + src_pt = mp.Vector3(-0.5*self.sx+self.dpml,0,0) + sources = [mp.Source(mp.GaussianSource(self.fcen,fwidth=self.df), + component=mp.Ez, + center=src_pt, + size=mp.Vector3(0,self.sy,0), + amp_func=pw_amp(k,src_pt))] + + sim = mp.Simulation(resolution=self.resolution, + cell_size=self.cell_size, + boundary_layers=self.abs_layers, + k_point=k, + default_material=self.glass, + sources=sources, + symmetries=symmetries, + kz_2d=kz_2d) + + refl_pt = mp.Vector3(-0.5*self.sx+self.dpml+0.5*self.dsub,0,0) + refl_flux = sim.add_mode_monitor(self.fcen, + 0, + 1, + mp.FluxRegion(center=refl_pt, + size=mp.Vector3(0,self.sy,0))) + + sim.run(until_after_sources=mp.stop_when_dft_decayed()) + + input_flux = mp.get_fluxes(refl_flux) + input_flux_data = sim.get_flux_data(refl_flux) + + sim.reset_meep() + + sim = mp.Simulation(resolution=self.resolution, + cell_size=self.cell_size, + boundary_layers=self.abs_layers, + geometry=self.geometry, + k_point=k, + sources=sources, + symmetries=symmetries, + kz_2d=kz_2d) + + refl_flux = sim.add_mode_monitor(self.fcen, + 0, + 1, + mp.FluxRegion(center=refl_pt, + size=mp.Vector3(0,self.sy,0))) + + sim.load_minus_flux_data(refl_flux,input_flux_data) + + tran_pt = mp.Vector3(0.5*self.sx-self.dpml-0.5*self.dpad,0,0) + tran_flux = sim.add_mode_monitor(self.fcen, + 0, + 1, + mp.FluxRegion(center=tran_pt, + size=mp.Vector3(0,self.sy,0))) + + sim.run(until_after_sources=mp.stop_when_dft_decayed()) + + # number of reflected orders + nm_r = (np.ceil((np.sqrt((self.fcen*self.ng)**2-k.z**2)-k.y)*self.gp) - + np.floor((-np.sqrt((self.fcen*self.ng)**2-k.z**2)-k.y)*self.gp)) + nm_r = int(nm_r/2) + + Rsum = 0 + for nm in range(nm_r): + for S_pol in [False, True]: + res = sim.get_eigenmode_coefficients(refl_flux, + mp.DiffractedPlanewave([0,nm,0], + mp.Vector3(1,0,0), + 1 if S_pol else 0, + 0 if S_pol else 1)) + r_coeffs = res.alpha + Rmode = abs(r_coeffs[0,0,1])**2/input_flux[0] + print("refl-order:, {}, {}, {}".format("s" if S_pol else "p",nm,Rmode)) + Rsum += Rmode if nm == 0 else 2*Rmode + + # number of transmitted orders + nm_t = (np.ceil((np.sqrt(self.fcen**2-k.z**2)-k.y)*self.gp) - + np.floor((-np.sqrt(self.fcen**2-k.z**2)-k.y)*self.gp)) + nm_t = int(nm_t/2) + + Tsum = 0 + for nm in range(nm_t): + for S_pol in [False, True]: + res = sim.get_eigenmode_coefficients(tran_flux, + mp.DiffractedPlanewave([0,nm,0], + mp.Vector3(1,0,0), + 1 if S_pol else 0, + 0 if S_pol else 1)) + t_coeffs = res.alpha + Tmode = abs(t_coeffs[0,0,0])**2/input_flux[0] + print("tran-order:, {}, {}, {}".format("s" if S_pol else "p",nm,Tmode)) + Tsum += Tmode if nm == 0 else 2*Tmode + + r_flux = mp.get_fluxes(refl_flux) + t_flux = mp.get_fluxes(tran_flux) + Rflux = -r_flux[0]/input_flux[0] + Tflux = t_flux[0]/input_flux[0] + + print("refl:, {}, {}".format(Rsum,Rflux)) + print("tran:, {}, {}".format(Tsum,Tflux)) + print("sum:, {}, {}".format(Rsum+Tsum,Rflux+Tflux)) + self.assertAlmostEqual(Rsum,Rflux,places=2) self.assertAlmostEqual(Tsum,Tflux,places=2) + self.assertAlmostEqual(Rsum+Tsum,1.00,places=2) - def test_binary_grating(self): - self.run_binary_grating_oblique(0) - self.run_binary_grating_oblique(10.7) if __name__ == '__main__': unittest.main() diff --git a/python/tests/test_special_kz.py b/python/tests/test_special_kz.py index ab119c8c3..8c0e051b2 100644 --- a/python/tests/test_special_kz.py +++ b/python/tests/test_special_kz.py @@ -1,26 +1,29 @@ import unittest +import parameterized import meep as mp import cmath import math from time import time + class TestSpecialKz(unittest.TestCase): def refl_planar(self, theta, kz_2d): resolution = 100 # pixels/um dpml = 1.0 - sx = 3+2*dpml + sx = 3.0 + 2*dpml sy = 1/resolution cell_size = mp.Vector3(sx,sy) pml_layers = [mp.PML(dpml,direction=mp.X)] - fcen = 1.0 # source wavelength = 1 um + fcen = 1.0 - k_point = mp.Vector3(z=math.sin(theta)).scale(fcen) + # plane of incidence is XZ + k_point = mp.Vector3(1,0,0).rotate(mp.Vector3(0,1,0),theta).scale(fcen) sources = [mp.Source(mp.GaussianSource(fcen,fwidth=0.2*fcen), - component=mp.Ez, + component=mp.Ez, # P-polarization center=mp.Vector3(-0.5*sx+dpml), size=mp.Vector3(y=sy))] @@ -63,6 +66,7 @@ def refl_planar(self, theta, kz_2d): Rmeep = -refl_flux[0]/empty_flux[0] return Rmeep + def test_special_kz(self): n1 = 1 n2 = 3.5 @@ -73,16 +77,18 @@ def test_special_kz(self): # compute Fresnel reflectance for P-polarization in medium n2 # for incident planewave in medium n1 at angle theta_in - Rfresnel = lambda theta_in: math.fabs((n1*math.cos(theta_out(theta_in))-n2*math.cos(theta_in))/(n1*math.cos(theta_out(theta_in))+n2*math.cos(theta_in)))**2 + Rfresnel = lambda theta_in: math.fabs((n1*math.cos(theta_out(theta_in))-n2*math.cos(theta_in)) + / (n1*math.cos(theta_out(theta_in))+n2*math.cos(theta_in)))**2 + # angle of incident planewave; clockwise (CW) about Y axis, 0 degrees along +X axis theta = math.radians(23) start = time() - Rmeep_complex = self.refl_planar(theta, 'complex') + Rmeep_complex = self.refl_planar(theta, "complex") t_complex = time() - start start = time() - Rmeep_real_imag = self.refl_planar(theta, 'real/imag') + Rmeep_real_imag = self.refl_planar(theta, "real/imag") t_real_imag = time() - start Rfres = Rfresnel(theta) @@ -90,14 +96,17 @@ def test_special_kz(self): self.assertAlmostEqual(Rmeep_complex,Rfres,places=2) self.assertAlmostEqual(Rmeep_real_imag,Rfres,places=2) - # the real/imag algorithm should be faster, but on CI machines performance is too variable for this to reliably hold + # the real/imag algorithm should be faster, but on CI machines performance is too variable + # for this to reliably hold # self.assertLess(t_real_imag,t_complex) - def eigsrc_kz(self, kz_2d): - print(kz_2d) + @parameterized.parameterized.expand([("complex",),("real/imag",)]) + def test_eigsrc_kz(self, kz_2d): resolution = 30 # pixels/um + cell_size = mp.Vector3(14,14) + pml_layers = [mp.PML(thickness=2)] geometry = [mp.Block(center=mp.Vector3(), @@ -124,8 +133,14 @@ def eigsrc_kz(self, kz_2d): k_point=mp.Vector3(z=kz), kz_2d=kz_2d) - tran = sim.add_flux(fsrc, 0, 1, mp.FluxRegion(center=mp.Vector3(x=5), size=mp.Vector3(y=14))) + tran = sim.add_flux(fsrc, + 0, + 1, + mp.FluxRegion(center=mp.Vector3(x=5), + size=mp.Vector3(y=14))) + sim.run(until_after_sources=50) + res = sim.get_eigenmode_coefficients(tran, [1,2], eig_parity=mp.EVEN_Y) @@ -146,10 +161,6 @@ def eigsrc_kz(self, kz_2d): self.assertAlmostEqual(ratio_ez.real,phase_diff.real,places=10) self.assertAlmostEqual(ratio_ez.imag,phase_diff.imag,places=10) - def test_eigsrc_kz(self): - self.eigsrc_kz("complex") - self.eigsrc_kz("real/imag") - if __name__ == '__main__': unittest.main() diff --git a/src/mpb.cpp b/src/mpb.cpp index d6afe8c66..63f701a44 100644 --- a/src/mpb.cpp +++ b/src/mpb.cpp @@ -292,6 +292,18 @@ static int nextpow2357(int n) { } } +void special_kz_phasefix(eigenmode_data *edata, bool phase_flip) { + size_t n = edata->n[0] * edata->n[1] * edata->n[2]; + complex *E = (complex *)edata->fft_data_E; + complex *H = (complex *)edata->fft_data_H; + complex im(0, phase_flip ? -1 : 1); + for (size_t i = 0; i < n; ++i) { + E[3*i + 2] *= im; // Ez + H[3*i + 0] *= im; // Hx + H[3*i + 1] *= im; // Hy + } +} + /****************************************************************/ /* call MPB to get the band_numth eigenmode at freq frequency. */ /* */ @@ -333,8 +345,8 @@ void *fields::get_eigenmode(double frequency, direction d, const volume where, c // special case: 2d cell in x and y with non-zero kz if ((eig_vol.dim == D3) && (float(v.in_direction(Z)) == float(1 / a)) && - (boundaries[High][Z] == Periodic && boundaries[Low][Z] == Periodic) && (kpoint.z() == 0) && - (real(k[Z]) != 0)) { + (boundaries[High][Z] == Periodic && boundaries[Low][Z] == Periodic) && + (kpoint.z() == 0) && (real(k[Z]) != 0)) { kpoint.set_direction(Z, real(k[Z])); empty_dim[2] = true; } @@ -609,11 +621,24 @@ void *fields::get_eigenmode(double frequency, direction d, const volume where, c k2sum += ktmp*ktmp; } } + if (((v.dim == D3) && (float(v.in_direction(Z)) == float(1 / a)) && + (boundaries[High][Z] == Periodic && boundaries[Low][Z] == Periodic) && (real(k[Z]) != 0)) || + ((v.dim == D2) && (real(k[Z]) != 0))) { + k2sum += k[Z]*k[Z]; + } // compute kperp (if it is non evanescent) OR // frequency from kperp^2 and sum of (kparallel+G)^2 { direction dd = eig_vol.normal_direction(); + if (eig_vol.dim == D3 && empty_dim[2]) { + volume eig_vol_2d(D2); + eig_vol_2d.set_direction_min(X, eig_vol.in_direction_min(X)); + eig_vol_2d.set_direction_min(Y, eig_vol.in_direction_min(Y)); + eig_vol_2d.set_direction_max(X, eig_vol.in_direction_max(X)); + eig_vol_2d.set_direction_max(Y, eig_vol.in_direction_max(Y)); + dd = eig_vol_2d.normal_direction(); + } if (match_frequency) { vec cen = eig_vol.center(); double nn = sqrt(real(get_eps(cen, frequency)) * real(get_mu(cen, frequency))); @@ -809,6 +834,8 @@ void fields::add_eigenmode_source(component c0, const src_time &src, direction d NULL, NULL, dp); finished_working(); + if (is_real && beta != 0) special_kz_phasefix(global_eigenmode_data, true /* phase_flip */); + if (global_eigenmode_data == NULL) meep::abort("MPB could not find the eigenmode"); /* add_volume_source amp_fun coordinates are relative to where.center(); @@ -923,6 +950,8 @@ void fields::get_eigenmode_coefficients(dft_flux flux, const volume &eig_vol, in continue; } + if (is_real && beta != 0) special_kz_phasefix((eigenmode_data *) mode_data, false /* phase_flip */); + double vg = get_group_velocity(mode_data); vec kfound = get_k(mode_data); if (vgrp) vgrp[nb * num_freqs + nf] = vg;