From 173675baf051229b750cd3291d93731fd97b6d64 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Thu, 14 Mar 2019 14:12:23 -0400 Subject: [PATCH 1/6] first attempt to hack in periodic green's functions --- src/meep.hpp | 6 +++++- src/near2far.cpp | 54 ++++++++++++++++++++++++++++++++++++++++++------ 2 files changed, 53 insertions(+), 7 deletions(-) diff --git a/src/meep.hpp b/src/meep.hpp index 9793c1185..95ec32185 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -1117,7 +1117,8 @@ class dft_near2far { /* fourier tranforms of tangential E and H field components in a medium with the given scalar eps and mu */ dft_near2far(dft_chunk *F, double fmin, double fmax, int Nf, double eps, double mu, - const volume &where_); + const volume &where_, const direction periodic_d_[2], + const int periodic_n_[2], const double periodic_k_[2], const double period_[2]); dft_near2far(const dft_near2far &f); /* return an array (Ex,Ey,Ez,Hx,Hy,Hz) x Nfreq of the far fields at x */ @@ -1157,6 +1158,9 @@ class dft_near2far { dft_chunk *F; double eps, mu; volume where; + direction periodic_d[2]; + int periodic_n[2]; + double periodic_k[2], period[2]; }; /* Class to compute local-density-of-states spectra: the power spectrum diff --git a/src/near2far.cpp b/src/near2far.cpp index e6f936f9a..52a5d25c2 100644 --- a/src/near2far.cpp +++ b/src/near2far.cpp @@ -30,16 +30,30 @@ using namespace std; namespace meep { dft_near2far::dft_near2far(dft_chunk *F_, double fmin, double fmax, int Nf, double eps_, double mu_, - const volume &where_) + const volume &where_, const direction periodic_d_[2], + const int periodic_n_[2], const double periodic_k_[2], const double period_[2]) : Nfreq(Nf), F(F_), eps(eps_), mu(mu_), where(where_) { if (Nf <= 1) fmin = fmax = (fmin + fmax) * 0.5; freq_min = fmin; dfreq = Nf <= 1 ? 0.0 : (fmax - fmin) / (Nf - 1); + for (int i = 0; i < 2; ++i) { + periodic_d[i] = periodic_d_[i]; + periodic_n[i] = periodic_n_[i]; + periodic_k[i] = periodic_k_[i]; + period[i] = period_[i]; + } } dft_near2far::dft_near2far(const dft_near2far &f) : freq_min(f.freq_min), dfreq(f.dfreq), Nfreq(f.Nfreq), F(f.F), eps(f.eps), mu(f.mu), - where(f.where) {} + where(f.where) { + for (int i = 0; i < 2; ++i) { + periodic_d[i] = f.periodic_d[i]; + periodic_n[i] = f.periodic_n[i]; + periodic_k[i] = f.periodic_k[i]; + period[i] = f.period[i]; + } +} void dft_near2far::remove() { while (F) { @@ -241,9 +255,21 @@ void dft_near2far::farfield_lowlevel(std::complex *EH, const vec &x) { x0 = f->S.transform(x0, f->sn) + rshift; for (int i = 0; i < Nfreq; ++i) { double freq = freq_min + i * dfreq; - green(EH6, x, freq, eps, mu, x0, c0, f->dft[Nfreq * idx_dft + i]); - for (int j = 0; j < 6; ++j) - EH[i * 6 + j] += EH6[j]; + vec xs(x0); + for (int i0 = -periodic_n[0]; i0 <= periodic_n[0]; ++i0) { + if (periodic_d[0] != NO_DIRECTION) + xs.set_direction(periodic_d[0], x0.in_direction(periodic_d[0]) + i0*period[0]); + double phase0 = i0*periodic_k[0]; + for (int i1 = -periodic_n[1]; i1 <= periodic_n[1]; ++i1) { + if (periodic_d[1] != NO_DIRECTION) + xs.set_direction(periodic_d[1], x0.in_direction(periodic_d[1]) + i1*period[1]); + double phase = phase0 + i1*periodic_k[1]; + std::complex cphase = std::polar(1.0, phase); + green(EH6, x, freq, eps, mu, x0, c0, f->dft[Nfreq * idx_dft + i]); + for (int j = 0; j < 6; ++j) + EH[i * 6 + j] += EH6[j] * cphase; + } + } } idx_dft++; } @@ -423,6 +449,10 @@ dft_near2far fields::add_dft_near2far(const volume_list *where, double freq_min, double eps = 0, mu = 0; volume everywhere = where->v; + direction periodic_d[2] = {NO_DIRECTION, NO_DIRECTION}; + int periodic_n[2] = {0, 0}; + double periodic_k[2] = {0, 0}, period[2] = {0, 0}; + for (const volume_list *w = where; w; w = w->next) { everywhere = everywhere | where->v; direction nd = component_direction(w->c); @@ -465,6 +495,17 @@ dft_near2far fields::add_dft_near2far(const volume_list *where, double freq_min, default: abort("invalid normal direction in dft_near2far!"); } + 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]))) { + periodic_d[i] = fd[i]; + periodic_n[i] = 10; // just hard-code a 10x10 supercell for now + period[i] = v.in_direction(fd[i]); + periodic_k[i] = 2*pi*real(k[fd[i]]) * period[i]; + } + } + for (int i = 0; i < 2; ++i) { /* E or H */ for (int j = 0; j < 2; ++j) { /* first or second component */ component c = direction_component(i == 0 ? Ex : Hx, fd[j]); @@ -480,7 +521,8 @@ dft_near2far fields::add_dft_near2far(const volume_list *where, double freq_min, } } - return dft_near2far(F, freq_min, freq_max, Nfreq, eps, mu, everywhere); + return dft_near2far(F, freq_min, freq_max, Nfreq, eps, mu, everywhere, + periodic_d, periodic_n, periodic_k, period); } } // namespace meep From 7492faff43d91e03e4c4c13432fdd53824220a38 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Thu, 21 Mar 2019 14:07:06 -0400 Subject: [PATCH 2/6] typo --- src/near2far.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/near2far.cpp b/src/near2far.cpp index 52a5d25c2..28c382d18 100644 --- a/src/near2far.cpp +++ b/src/near2far.cpp @@ -265,7 +265,7 @@ void dft_near2far::farfield_lowlevel(std::complex *EH, const vec &x) { xs.set_direction(periodic_d[1], x0.in_direction(periodic_d[1]) + i1*period[1]); double phase = phase0 + i1*periodic_k[1]; std::complex cphase = std::polar(1.0, phase); - green(EH6, x, freq, eps, mu, x0, c0, f->dft[Nfreq * idx_dft + i]); + green(EH6, x, freq, eps, mu, xs, c0, f->dft[Nfreq * idx_dft + i]); for (int j = 0; j < 6; ++j) EH[i * 6 + j] += EH6[j] * cphase; } From e719275c9c2774c538bfe1e7e9aac1a85827ad5c Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Thu, 21 Mar 2019 14:25:05 -0400 Subject: [PATCH 3/6] make nperiods an option for near2far --- python/simulation.py | 19 ++++++++++++------- src/meep.hpp | 2 +- src/near2far.cpp | 7 +++---- 3 files changed, 16 insertions(+), 12 deletions(-) diff --git a/python/simulation.py b/python/simulation.py index f7b947104..a3540b713 100644 --- a/python/simulation.py +++ b/python/simulation.py @@ -321,7 +321,8 @@ class DftNear2Far(DftObj): def __init__(self, func, args): super(DftNear2Far, self).__init__(func, args) self.nfreqs = args[2] - self.regions = args[3] + self.nperiods = args[3] + self.regions = args[4] self.num_components = 4 @property @@ -1371,15 +1372,19 @@ def get_dft_data(self, dft_chunk): mp._get_dft_data(dft_chunk, arr) return arr - def add_near2far(self, fcen, df, nfreq, *near2fars): - n2f = DftNear2Far(self._add_near2far, [fcen, df, nfreq, near2fars]) + def add_near2far(self, fcen, df, nfreq, nperiods=1, *near2fars): + # hack around python 2's inability to put keyword args after varargs: + if not isinstance(nperiods, int): + near2fars = (nperiods,) + near2fars + nperiods = 1 + n2f = DftNear2Far(self._add_near2far, [fcen, df, nfreq, nperiods, near2fars]) self.dft_objects.append(n2f) return n2f - def _add_near2far(self, fcen, df, nfreq, near2fars): + def _add_near2far(self, fcen, df, nfreq, nperiods, near2fars): if self.fields is None: self.init_sim() - return self._add_fluxish_stuff(self.fields.add_dft_near2far, fcen, df, nfreq, near2fars) + return self._add_fluxish_stuff(self.fields.add_dft_near2far, fcen, df, nfreq, near2fars, nperiods) def add_energy(self, fcen, df, nfreq, *energys): en = DftEnergy(self._add_energy, [fcen, df, nfreq, energys]) @@ -1634,7 +1639,7 @@ def solve_cw(self, tol=1e-8, maxiters=10000, L=2): self._evaluate_dft_objects() return self.fields.solve_cw(tol, maxiters, L) - def _add_fluxish_stuff(self, add_dft_stuff, fcen, df, nfreq, stufflist): + def _add_fluxish_stuff(self, add_dft_stuff, fcen, df, nfreq, stufflist, *args): vol_list = None for s in stufflist: @@ -1647,7 +1652,7 @@ def _add_fluxish_stuff(self, add_dft_stuff, fcen, df, nfreq, stufflist): is_cylindrical=self.is_cylindrical).swigobj vol_list = mp.make_volume_list(v2, c, s.weight, vol_list) - stuff = add_dft_stuff(vol_list, fcen - df / 2, fcen + df / 2, nfreq) + stuff = add_dft_stuff(vol_list, fcen - df / 2, fcen + df / 2, nfreq, *args) vol_list.__swig_destroy__(vol_list) return stuff diff --git a/src/meep.hpp b/src/meep.hpp index 95ec32185..82d7a0626 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -1733,7 +1733,7 @@ class fields { // near2far.cpp dft_near2far add_dft_near2far(const volume_list *where, double freq_min, double freq_max, - int Nfreq); + int Nfreq, int Nperiods=1); // monitor.cpp double get_chi1inv(component, direction, const vec &loc) const; double get_inveps(component c, direction d, const vec &loc) const { diff --git a/src/near2far.cpp b/src/near2far.cpp index 28c382d18..1ce20f894 100644 --- a/src/near2far.cpp +++ b/src/near2far.cpp @@ -444,7 +444,7 @@ double *dft_near2far::flux(direction df, const volume &where, double resolution) static double approxeq(double a, double b) { return fabs(a - b) < 0.5e-11 * (fabs(a) + fabs(b)); } dft_near2far fields::add_dft_near2far(const volume_list *where, double freq_min, double freq_max, - int Nfreq) { + int Nfreq, int Nperiods) { dft_chunk *F = 0; /* E and H chunks*/ double eps = 0, mu = 0; volume everywhere = where->v; @@ -497,10 +497,9 @@ dft_near2far fields::add_dft_near2far(const volume_list *where, double freq_min, 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]))) { + boundaries[High][fd[i]] == Periodic && boundaries[Low][fd[i]] == Periodic) { periodic_d[i] = fd[i]; - periodic_n[i] = 10; // just hard-code a 10x10 supercell for now + periodic_n[i] = Nperiods; period[i] = v.in_direction(fd[i]); periodic_k[i] = 2*pi*real(k[fd[i]]) * period[i]; } From 8c75814c7eb745efe92231ff85188046ffa65eaf Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Thu, 21 Mar 2019 14:29:21 -0400 Subject: [PATCH 4/6] only do lattice sum in non-empty dirs --- src/near2far.cpp | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/src/near2far.cpp b/src/near2far.cpp index 1ce20f894..ba0547168 100644 --- a/src/near2far.cpp +++ b/src/near2far.cpp @@ -495,13 +495,16 @@ dft_near2far fields::add_dft_near2far(const volume_list *where, double freq_min, default: abort("invalid normal direction in dft_near2far!"); } - for (int i = 0; i < 2; ++i) { - if (has_direction(v.dim, fd[i]) && - boundaries[High][fd[i]] == Periodic && boundaries[Low][fd[i]] == Periodic) { - periodic_d[i] = fd[i]; - periodic_n[i] = Nperiods; - period[i] = v.in_direction(fd[i]); - periodic_k[i] = 2*pi*real(k[fd[i]]) * period[i]; + if (Nperiods > 1) { + 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])) >= 0.5*float(v.in_direction(fd[i]))) { + periodic_d[i] = fd[i]; + periodic_n[i] = Nperiods; + period[i] = v.in_direction(fd[i]); + periodic_k[i] = 2*pi*real(k[fd[i]]) * period[i]; + } } } From 93c3b3ad583269f99788abd2ee6fd120b9a01bdc Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Thu, 21 Mar 2019 14:33:12 -0400 Subject: [PATCH 5/6] simplify kw handling in add_near2far --- python/simulation.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/python/simulation.py b/python/simulation.py index a3540b713..a625dc117 100644 --- a/python/simulation.py +++ b/python/simulation.py @@ -1372,11 +1372,8 @@ def get_dft_data(self, dft_chunk): mp._get_dft_data(dft_chunk, arr) return arr - def add_near2far(self, fcen, df, nfreq, nperiods=1, *near2fars): - # hack around python 2's inability to put keyword args after varargs: - if not isinstance(nperiods, int): - near2fars = (nperiods,) + near2fars - nperiods = 1 + def add_near2far(self, fcen, df, nfreq, *near2fars, **kwargs): + nperiods = kwargs.get('nperiods', 1) n2f = DftNear2Far(self._add_near2far, [fcen, df, nfreq, nperiods, near2fars]) self.dft_objects.append(n2f) return n2f From 75cf9eaa88a4dec1f57a21cd17ac737d0e2c7066 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Mon, 25 Mar 2019 17:10:48 -0400 Subject: [PATCH 6/6] fix symmetry again --- src/near2far.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/near2far.cpp b/src/near2far.cpp index ba0547168..576e3d70c 100644 --- a/src/near2far.cpp +++ b/src/near2far.cpp @@ -497,12 +497,13 @@ dft_near2far fields::add_dft_near2far(const volume_list *where, double freq_min, if (Nperiods > 1) { for (int i = 0; i < 2; ++i) { + double user_width = user_volume.num_direction(fd[i]) / a; if (has_direction(v.dim, fd[i]) && boundaries[High][fd[i]] == Periodic && boundaries[Low][fd[i]] == Periodic && - float(w->v.in_direction(fd[i])) >= 0.5*float(v.in_direction(fd[i]))) { + float(w->v.in_direction(fd[i])) >= float(user_width)) { periodic_d[i] = fd[i]; periodic_n[i] = Nperiods; - period[i] = v.in_direction(fd[i]); + period[i] = user_width; periodic_k[i] = 2*pi*real(k[fd[i]]) * period[i]; } }