diff --git a/.gitignore b/.gitignore index 141f6092d..0b1fd93e5 100644 --- a/.gitignore +++ b/.gitignore @@ -73,6 +73,7 @@ tests/h5test tests/harmonics tests/integrate tests/known_results +tests/near2far tests/one_dimensional tests/physical tests/pml diff --git a/libctl/meep-ctl-swig.hpp b/libctl/meep-ctl-swig.hpp index 9b093b52c..6595e00c8 100644 --- a/libctl/meep-ctl-swig.hpp +++ b/libctl/meep-ctl-swig.hpp @@ -31,6 +31,7 @@ ctlio::number_list dft_force_force(meep::dft_force *f); ctlio::number_list dft_ldos_ldos(meep::dft_ldos *f); ctlio::cnumber_list dft_ldos_F(meep::dft_ldos *f); ctlio::cnumber_list dft_ldos_J(meep::dft_ldos *f); +ctlio::cnumber_list dft_near2far_farfield(meep::dft_near2far *f, const meep::vec &x); ctlio::cnumber_list make_casimir_g(double T, double dt, double sigma, meep::field_type ft, std::complex (*eps_func)(std::complex omega) = 0, diff --git a/libctl/meep.cpp b/libctl/meep.cpp index 98e857b46..0fe93251c 100644 --- a/libctl/meep.cpp +++ b/libctl/meep.cpp @@ -118,6 +118,14 @@ ctlio::cnumber_list dft_ldos_J(dft_ldos *f) return res; } +ctlio::cnumber_list dft_near2far_farfield(dft_near2far *f, const vec &x) +{ + ctlio::cnumber_list res; + res.num_items = f->Nfreq * 6; + res.items = (cnumber *) f->farfield(x); + return res; +} + /***************************************************************************/ ctlio::cnumber_list make_casimir_g(double T, double dt, double sigma, meep::field_type ft, diff --git a/libctl/meep.scm.in b/libctl/meep.scm.in index 84f1a9996..1b8f36728 100644 --- a/libctl/meep.scm.in +++ b/libctl/meep.scm.in @@ -625,6 +625,51 @@ (load-force fname force) (meep-dft-force-scale-dfts force -1.0)) +; **************************************************************** +; Near-to-far-field transformations (again similar to dft-foobar) + +(define-class near2far-region no-parent + (define-property center no-default 'vector3) + (define-property size (vector3 0 0 0) 'vector3) + (define-property direction AUTOMATIC 'integer) + (define-property weight 1.0 'cnumber)) + +(define (fields-add-near2far fields fcen df nfreq . near2fars) + (fields-add-fluxish-stuff meep-fields-add-dft-near2far + fields fcen df nfreq near2fars)) + +(define (add-near2far fcen df nfreq . near2fars) + (if (null? fields) (init-fields)) + (apply fields-add-near2far (append (list fields fcen df nfreq) near2fars))) + +(define (scale-near2far-fields s f) + (meep-dft-near2far-scale-dfts f s)) + +(define (get-near2far-freqs f) + (arith-sequence + (meep-dft-near2far-freq-min-get f) + (meep-dft-near2far-dfreq-get f) + (meep-dft-near2far-Nfreq-get f))) + +(define (get-farfield f x) + (dft-near2far-farfield f x)) + +(define (output-farfields near2far fname where resolution) + (meep-dft-near2far-save-farfields near2far fname (get-filename-prefix) + where resolution)) + +(define (load-near2far fname near2far) + (if (null? fields) (init-fields)) + (meep-dft-near2far-load-hdf5 near2far fields fname "" (get-filename-prefix))) + +(define (save-near2far fname near2far) + (if (null? fields) (init-fields)) + (meep-dft-near2far-save-hdf5 near2far fields fname "" (get-filename-prefix))) + +(define (load-minus-near2far fname near2far) + (load-near2far fname near2far) + (meep-dft-near2far-scale-dfts near2far -1.0)) + ; **************************************************************** ; Generic step functions: these are functions which are called ; (potentially) at every time step. They can either be a thunk diff --git a/src/Makefile.am b/src/Makefile.am index 29c11db9d..e77a35231 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -13,7 +13,7 @@ libmeep@MEEP_SUFFIX@_la_SOURCES = anisotropic_averaging.cpp bands.cpp \ boundaries.cpp bicgstab.cpp casimir.cpp control_c.cpp cw_fields.cpp \ dft.cpp dft_ldos.cpp energy_and_flux.cpp fields.cpp loop_in_chunks.cpp \ grace.cpp h5fields.cpp h5file.cpp initialize.cpp integrate.cpp \ -integrate2.cpp monitor.cpp mympi.cpp multilevel-atom.cpp \ +integrate2.cpp monitor.cpp mympi.cpp multilevel-atom.cpp near2far.cpp \ output_directory.cpp random.cpp sources.cpp step.cpp step_db.cpp \ stress.cpp structure.cpp susceptibility.cpp time.cpp update_eh.cpp \ mpb.cpp update_pols.cpp vec.cpp step_generic.cpp $(HDRS) \ diff --git a/src/cw_fields.cpp b/src/cw_fields.cpp index 925d696cb..b05c8d111 100644 --- a/src/cw_fields.cpp +++ b/src/cw_fields.cpp @@ -116,6 +116,7 @@ bool fields::solve_cw(double tol, int maxiters, complex frequency, int L) { if (is_real) abort("solve_cw is incompatible with use_real_fields()"); if (L < 1) abort("solve_cw called with L = %d < 1", L); + int tsave = t; // save time (gets incremented by iterations) set_solve_cw_omega(2*pi*frequency); @@ -145,8 +146,6 @@ bool fields::solve_cw(double tol, int maxiters, complex frequency, complex *x = reinterpret_cast*>(work + nwork); complex *b = reinterpret_cast*>(work + nwork + N); - int tsave = t; // save time (gets incremented by iterations) - fields_to_array(*this, x); // initial guess = initial fields // get J amplitudes from current time step diff --git a/src/dft.cpp b/src/dft.cpp index ba0ccf580..bcb253407 100644 --- a/src/dft.cpp +++ b/src/dft.cpp @@ -30,6 +30,7 @@ struct dft_chunk_data { // for passing to field::loop_in_chunks as void* double omega_min, domega; int Nomega; component c; + int vc; complex weight, extra_weight; bool include_dV_and_interp_weights; bool sqrt_dV_and_interp_weights; @@ -44,6 +45,7 @@ dft_chunk::dft_chunk(fields_chunk *fc_, complex scale_, component c_, bool use_centered_grid, + ivec shift_, const symmetry &S_, int sn_, int vc_, const void *data_) { dft_chunk_data *data = (dft_chunk_data *) data_; if (!fc_->f[c_][0]) @@ -84,6 +86,10 @@ dft_chunk::dft_chunk(fields_chunk *fc_, else avg1 = avg2 = 0; + shift = shift_; + S = S_; sn = sn_; + vc = vc_; + omega_min = data->omega_min; domega = data->domega; Nomega = data->Nomega; @@ -140,7 +146,7 @@ static void add_dft_chunkloop(fields_chunk *fc, int ichunk, component cgrid, void *chunkloop_data) { dft_chunk_data *data = (dft_chunk_data *) chunkloop_data; - (void) shift; (void) ichunk; // unused + (void) ichunk; // unused component c = S.transform(data->c, -sn); if (c >= NUM_FIELD_COMPONENTS || !fc->f[c][0]) @@ -150,6 +156,7 @@ static void add_dft_chunkloop(fields_chunk *fc, int ichunk, component cgrid, data->extra_weight, shift_phase * S.phase_shift(c, sn), c, cgrid == Centered, + shift, S, sn, data->vc, chunkloop_data); } @@ -159,12 +166,13 @@ dft_chunk *fields::add_dft(component c, const volume &where, complex weight, dft_chunk *chunk_next, bool sqrt_dV_and_interp_weights, complex extra_weight, - bool use_centered_grid) { + bool use_centered_grid, int vc) { if (coordinate_mismatch(gv.dim, c)) return NULL; dft_chunk_data data; data.c = c; + data.vc = vc; if (Nfreq <= 1) freq_min = freq_max = (freq_min + freq_max) * 0.5; data.omega_min = freq_min * 2*pi; data.domega = Nfreq <= 1 ? 0.0 : diff --git a/src/meep.hpp b/src/meep.hpp index b641a5b97..2e31c1482 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -827,6 +827,7 @@ class dft_chunk { std::complex extra_weight_, component c_, bool use_centered_grid, + ivec shift_, const symmetry &S_, int sn_, int vc, const void *data_); ~dft_chunk(); @@ -857,7 +858,6 @@ class dft_chunk { it is used in computations involving dft[...], it needs to be public. */ std::complex extra_weight; -private: // parameters passed from field_integrate: fields_chunk *fc; ivec is, ie; @@ -865,11 +865,15 @@ class dft_chunk { double dV0, dV1; bool sqrt_dV_and_interp_weights; std::complex scale; // scale factor * phase from shift and symmetry + ivec shift; + symmetry S; int sn; // cache of exp(iwt) * scale, of length Nomega std::complex *dft_phase; int avg1, avg2; // index offsets for average to get epsilon grid + + int vc; // component descriptor from the original volume }; void save_dft_hdf5(dft_chunk *dft_chunks, component c, h5file *file, @@ -939,6 +943,48 @@ class dft_force { dft_chunk *offdiag1, *offdiag2, *diag; }; +// near2far.cpp (normally created with fields::add_dft_near2far) +class dft_near2far { +public: + /* 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); + dft_near2far(const dft_near2far &f); + + /* return an array (Ex,Ey,Ez,Hx,Hy,Hz) x Nfreq of the far fields at x */ + std::complex *farfield(const vec &x); + + /* like farfield, but requires F to be Nfreq*6 preallocated array, and + does *not* perform the reduction over processes...an MPI allreduce + summation by the caller is required to get the final result ... used + by other output routine to efficiently get far field on a grid of pts */ + void farfield_lowlevel(std::complex *F, const vec &x); + + /* output far fields on a grid to an HDF5 file */ + void save_farfields(const char *fname, const char *prefix, + const volume &where, double resolution); + + void save_hdf5(h5file *file, const char *dprefix = 0); + void load_hdf5(h5file *file, const char *dprefix = 0); + + void operator-=(const dft_near2far &fl); + + void save_hdf5(fields &f, const char *fname, const char *dprefix = 0, + const char *prefix = 0); + void load_hdf5(fields &f, const char *fname, const char *dprefix = 0, + const char *prefix = 0); + + void scale_dfts(std::complex scale); + + void remove(); + + double freq_min, dfreq; + int Nfreq; + dft_chunk *F; + double eps, mu; +}; + /* Class to compute local-density-of-states spectra: the power spectrum P(omega) of the work done by the sources. Specialized to handle only the case where all sources have the same time dependence, which greatly @@ -1348,7 +1394,7 @@ class fields { std::complex weight = 1.0, dft_chunk *chunk_next = 0, bool sqrt_dV_and_interp_weights = false, std::complex extra_weight = 1.0, - bool use_centered_grid = true); + bool use_centered_grid = true, int vc = 0); dft_chunk *add_dft_pt(component c, const vec &where, double freq_min, double freq_max, int Nfreq); dft_chunk *add_dft(const volume_list *where, @@ -1368,6 +1414,9 @@ class fields { dft_force add_dft_force(const volume_list *where, double freq_min, double freq_max, int Nfreq); + // near2far.cpp + dft_near2far add_dft_near2far(const volume_list *where, + double freq_min, double freq_max, int Nfreq); // monitor.cpp double get_chi1inv(component, direction, const vec &loc) const; double get_inveps(component c, direction d, const vec &loc) const { @@ -1570,6 +1619,16 @@ int random_int(int a, int b); // uniform random in [a,b) // Bessel function (in initialize.cpp) double BesselJ(int m, double kr); +// analytical Green's functions (in near2far.cpp); upon return, +// EH[0..5] are set to the Ex,Ey,Ez,Hx,Hy,Hz field components at x +// from a c0 source of amplitude f0 at x0. +void green2d(std::complex *EH, const vec &x, + double freq, double eps, double mu, + const vec &x0, component c0, std::complex f0); +void green3d(std::complex *EH, const vec &x, + double freq, double eps, double mu, + const vec &x0, component c0, std::complex f0); + } /* namespace meep */ #endif /* MEEP_H */ diff --git a/src/meep/mympi.hpp b/src/meep/mympi.hpp index 432e94361..c7e7d4fe8 100644 --- a/src/meep/mympi.hpp +++ b/src/meep/mympi.hpp @@ -65,9 +65,11 @@ int max_to_all(int); double sum_to_master(double); // Only returns the correct value to proc 0. double sum_to_all(double); void sum_to_all(const double *in, double *out, int size); +void sum_to_master(const double *in, double *out, int size); void sum_to_all(const float *in, double *out, int size); void sum_to_all(const std::complex *in, std::complex *out, int size); void sum_to_all(const std::complex *in, std::complex *out, int size); +void sum_to_master(const std::complex *in, std::complex *out, int size); long double sum_to_all(long double); std::complex sum_to_all(std::complex in); std::complex sum_to_all(std::complex in); diff --git a/src/meep/vec.hpp b/src/meep/vec.hpp index 03a8431aa..d0b4c10b5 100644 --- a/src/meep/vec.hpp +++ b/src/meep/vec.hpp @@ -918,6 +918,14 @@ class volume_list { public: volume_list(const volume &v, int c, std::complex weight = 1.0, volume_list *next = 0) : v(v), c(c), weight(weight), next(next) {} ~volume_list() { delete next; } + volume_list(const volume_list *vl) : v(vl->v), c(vl->c), weight(vl->weight), next(0) { + volume_list *p = vl->next, *q = this; + while (p) { + q->next = new volume_list(*p); + q = q->next; + p = p->next; + } + } volume v; int c; // component or derived component associated with v (e.g. for flux) diff --git a/src/mympi.cpp b/src/mympi.cpp index 385e934a1..a4200200c 100644 --- a/src/mympi.cpp +++ b/src/mympi.cpp @@ -293,6 +293,14 @@ void sum_to_all(const double *in, double *out, int size) { #endif } +void sum_to_master(const double *in, double *out, int size) { +#ifdef HAVE_MPI + MPI_Reduce((void*) in, out, size, MPI_DOUBLE,MPI_SUM,0,mycomm); +#else + memcpy(out, in, sizeof(double) * size); +#endif +} + void sum_to_all(const float *in, double *out, int size) { double *in2 = new double[size]; for (int i = 0; i < size; ++i) in2[i] = in[i]; @@ -308,6 +316,10 @@ void sum_to_all(const complex *in, complex *out, int size) { sum_to_all((const float*) in, (double*) out, 2*size); } +void sum_to_master(const complex *in, complex *out, int size) { + sum_to_master((const double*) in, (double*) out, 2*size); +} + long double sum_to_all(long double in) { long double out = in; #ifdef HAVE_MPI diff --git a/src/near2far.cpp b/src/near2far.cpp new file mode 100644 index 000000000..c5cf7b628 --- /dev/null +++ b/src/near2far.cpp @@ -0,0 +1,420 @@ +/* Copyright (C) 2005-2014 Massachusetts Institute of Technology. + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + */ + +/* Near-to-far field transformation: compute DFT of tangential fields on + a "near" surface, and use these (via the equivalence principle) to + compute the fields on a "far" surface via the homogeneous-medium Green's + function in 2d or 3d. */ + +#include +#include +#include "config.h" + +using namespace std; + +namespace meep { + +dft_near2far::dft_near2far(dft_chunk *F_, + double fmin, double fmax, int Nf, + double eps_, double mu_) +{ + if (Nf <= 1) fmin = fmax = (fmin + fmax) * 0.5; + freq_min = fmin; + Nfreq = Nf; + dfreq = Nf <= 1 ? 0.0 : (fmax - fmin) / (Nf - 1); + F = F_; + eps = eps_; mu = mu_; +} + +dft_near2far::dft_near2far(const dft_near2far &f) { + freq_min = f.freq_min; Nfreq = f.Nfreq; dfreq = f.dfreq; + F = f.F; + eps = f.eps; mu = f.mu; +} + +void dft_near2far::remove() +{ + while (F) { + dft_chunk *nxt = F->next_in_dft; + delete F; + F = nxt; + } +} + +void dft_near2far::operator-=(const dft_near2far &st) { + if (F && st.F) *F -= *st.F; +} + +void dft_near2far::save_hdf5(h5file *file, const char *dprefix) { + save_dft_hdf5(F, "F", file, dprefix); +} + +void dft_near2far::load_hdf5(h5file *file, const char *dprefix) { + load_dft_hdf5(F, "F", file, dprefix); +} + +void dft_near2far::save_hdf5(fields &f, const char *fname, const char *dprefix, + const char *prefix) { + h5file *ff = f.open_h5file(fname, h5file::WRITE, prefix); + save_hdf5(ff, dprefix); + delete ff; +} + +void dft_near2far::load_hdf5(fields &f, const char *fname, const char *dprefix, + const char *prefix) { + h5file *ff = f.open_h5file(fname, h5file::READONLY, prefix); + load_hdf5(ff, dprefix); + delete ff; +} + +void dft_near2far::scale_dfts(complex scale) { + if (F) F->scale_dft(scale); +} + +typedef void (*greenfunc)(std::complex *EH, const vec &x, + double freq, double eps, double mu, + const vec &x0, component c0, std::complex); + +/* Given the field f0 correponding to current-source component c0 at + x0, compute the E/H fields EH[6] (6 components) at x for a frequency + freq in the homogeneous 3d medium eps and mu. + + Adapted from code by M. T. Homer Reid in his SCUFF-EM package + (file scuff-em/src/libs/libIncField/PointSource.cc), which is GPL v2+. */ +void green3d(std::complex *EH, const vec &x, + double freq, double eps, double mu, + const vec &x0, component c0, std::complex f0) +{ + vec rhat = x - x0; + double r = abs(rhat); + rhat = rhat / r; + + if (rhat.dim != D3) abort("wrong dimensionality in green3d"); + + double n = sqrt(eps*mu); + double k = 2*pi*freq*n; + std::complex ikr = std::complex(0.0, k*r); + double ikr2 = -(k*r)*(k*r); + /* note that SCUFF-EM computes the fields from the dipole moment p, + whereas we need it from the current J = -i*omega*p, so our result + is divided by -i*omega compared to SCUFF */ + std::complex expfac = f0 * polar(k*n/(4*pi*r), k*r + pi*0.5); + double Z = sqrt(mu/eps); + + vec p = zero_vec(rhat.dim); + p.set_direction(component_direction(c0), 1); + double pdotrhat = p & rhat; + vec rhatcrossp = vec(rhat.y() * p.z() - + rhat.z() * p.y(), + rhat.z() * p.x() - + rhat.x() * p.z(), + rhat.x() * p.y() - + rhat.y() * p.x()); + + /* compute the various scalar quantities in the point source formulae */ + std::complex term1 = 1.0 - 1.0/ikr + 1.0/ikr2; + std::complex term2 = (-1.0 + 3.0/ikr - 3.0/ikr2) * pdotrhat; + std::complex term3 = (1.0 - 1.0/ikr); + + /* now assemble everything based on source type */ + if (is_electric(c0)) { + expfac /= eps; + + EH[0] = expfac * (term1*p.x() + term2*rhat.x()); + EH[1] = expfac * (term1*p.y() + term2*rhat.y()); + EH[2] = expfac * (term1*p.z() + term2*rhat.z()); + + EH[3] = expfac*term3*rhatcrossp.x() / Z; + EH[4] = expfac*term3*rhatcrossp.y() / Z; + EH[5] = expfac*term3*rhatcrossp.z() / Z; + } + else if (is_magnetic(c0)) { + expfac /= mu; + + EH[0] = -expfac*term3*rhatcrossp.x() * Z; + EH[1] = -expfac*term3*rhatcrossp.y() * Z; + EH[2] = -expfac*term3*rhatcrossp.z() * Z; + + EH[3] = expfac * (term1*p.x() + term2*rhat.x()); + EH[4] = expfac * (term1*p.y() + term2*rhat.y()); + EH[5] = expfac * (term1*p.z() + term2*rhat.z()); + } + else + abort("unrecognized source type"); +} + +#ifdef HAVE_LIBGSL +# include +// hankel function J + iY +static std::complex hankel(int n, double x) { + return std::complex(gsl_sf_bessel_Jn(n, x), + gsl_sf_bessel_Yn(n, x)); +} +#else /* !HAVE_LIBGSL */ +static std::complex hankel(int n, double x) { + (void) n; (void) x; // unused + abort("GNU GSL library is required for Hankel functions"); +} +#endif /* !HAVE_LIBGSL */ + +/* like green3d, but 2d Green's functions */ +void green2d(std::complex *EH, const vec &x, + double freq, double eps, double mu, + const vec &x0, component c0, std::complex f0) +{ + vec rhat = x - x0; + double r = abs(rhat); + rhat = rhat / r; + + if (rhat.dim != D2) abort("wrong dimensionality in green2d"); + + double omega = 2*pi*freq; + double k = omega*sqrt(eps*mu); + std::complex ik = std::complex(0.0, k); + double kr = k*r; + double Z = sqrt(mu/eps); + std::complex H0 = hankel(0, kr) * f0; + std::complex H1 = hankel(1, kr) * f0; + std::complex ikH1 = 0.25 * ik * H1; + + if (component_direction(c0) == meep::Z) { + if (is_electric(c0)) { // Ez source + EH[0] = EH[1] = 0.0; + EH[2] = (-0.25*omega*mu) * H0; + + EH[3] = -rhat.y() * ikH1; + EH[4] = rhat.x() * ikH1; + EH[5] = 0.0; + } + else /* (is_magnetic(c0)) */ { // Hz source + EH[0] = rhat.y() * ikH1; + EH[1] = -rhat.x() * ikH1; + EH[2] = 0.0; + + EH[3] = EH[4] = 0.0; + EH[5] = (-0.25*omega*eps) * H0; + } + } + else { /* in-plane source */ + std::complex H2 = hankel(2, kr) * f0; + + vec p = zero_vec(rhat.dim); + p.set_direction(component_direction(c0), 1); + + double pdotrhat = p & rhat; + double rhatcrossp = rhat.x() * p.y() - rhat.y() * p.x(); + + if (is_electric(c0)) { // Exy source + EH[0] = -(rhat.x() * (pdotrhat/r * 0.25*Z)) * H1 + + (rhat.y() * (rhatcrossp * omega*mu * 0.125)) * (H0 - H2); + EH[1] = -(rhat.y() * (pdotrhat/r * 0.25*Z)) * H1 - + (rhat.x() * (rhatcrossp * omega*mu * 0.125)) * (H0 - H2); + EH[2] = 0.0; + + EH[3] = EH[4] = 0.0; + EH[5] = -rhatcrossp * ikH1; + } + else /* (is_magnetic(c0)) */ { // Hxy source + EH[0] = EH[1] = 0.0; + EH[2] = rhatcrossp * ikH1; + + EH[3] = -(rhat.x() * (pdotrhat/r * 0.25/Z)) * H1 + + (rhat.y() * (rhatcrossp * omega*eps * 0.125)) * (H0 - H2); + EH[4] = -(rhat.y() * (pdotrhat/r * 0.25/Z)) * H1 - + (rhat.x() * (rhatcrossp * omega*eps * 0.125)) * (H0 - H2); + EH[5] = 0.0; + } + } +} + +void dft_near2far::farfield_lowlevel(std::complex *EH, const vec &x) +{ + if (x.dim != D3 && x.dim != D2) + abort("only 2d or 3d far-field computation is supported"); + greenfunc green = x.dim == D2 ? green2d : green3d; + + std::complex EH6[6]; + for (int i = 0; i < 6 * Nfreq; ++i) + EH[i] = 0.0; + + for (dft_chunk *f = F; f; f = f->next_in_dft) { + assert(Nfreq == f->Nomega); + + component c0 = component(f->vc); /* equivalent source component */ + + vec rshift(f->shift * (0.5*f->fc->gv.inva)); + int idx_dft = 0; + LOOP_OVER_IVECS(f->fc->gv, f->is, f->ie, idx) { + IVEC_LOOP_LOC(f->fc->gv, x0); + 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]; + } + idx_dft++; + } + } +} + +std::complex *dft_near2far::farfield(const vec &x) { + std::complex *EH, *EH_local; + EH_local = new std::complex[6 * Nfreq]; + farfield_lowlevel(EH_local, x); + EH = new std::complex[6 * Nfreq]; + sum_to_all(EH_local, EH, 6 * Nfreq); + delete[] EH_local; + return EH; +} + +void dft_near2far::save_farfields(const char *fname, const char *prefix, + const volume &where, double resolution) { + /* compute output grid size etc. */ + int dims[4] = {1,1,1,1}; + double dx[3] = {0,0,0}; + direction dirs[3] = {X,Y,Z}; + int rank = 0, N = 1; + LOOP_OVER_DIRECTIONS(where.dim, d) { + dims[rank] = int(floor(where.in_direction(d) * resolution)); + if (dims[rank] <= 1) { + dims[rank] = 1; + dx[rank] = 0; + } + else + dx[rank] = where.in_direction(d) / (dims[rank] - 1); + N *= dims[rank]; + dirs[rank++] = d; + } + + if (N * Nfreq < 1) return; /* nothing to output */ + + /* 6 x 2 x N x Nfreq array of fields in row-major order */ + realnum *EH = new realnum[6*2*N*Nfreq]; + realnum *EH_ = new realnum[6*2*N*Nfreq]; // temp array for sum_to_master + + /* fields for farfield_lowlevel for a single output point x */ + std::complex *EH1 = new std::complex[6*Nfreq]; + + vec x(where.dim); + for (int i0 = 0; i0 < dims[0]; ++i0) { + x.set_direction(dirs[0], where.in_direction_min(dirs[0]) + i0*dx[0]); + for (int i1 = 0; i1 < dims[1]; ++i1) { + x.set_direction(dirs[1], + where.in_direction_min(dirs[1]) + i1*dx[1]); + for (int i2 = 0; i2 < dims[2]; ++i2) { + x.set_direction(dirs[2], + where.in_direction_min(dirs[2]) + i2*dx[2]); + farfield_lowlevel(EH1, x); + int idx = (i0 * dims[1] + i1) * dims[2] + i2; + for (int i = 0; i < Nfreq; ++i) + for (int k = 0; k < 6; ++k) { + EH_[((k * 2 + 0) * N + idx) * Nfreq + i] = + real(EH1[i * 6 + k]); + EH_[((k * 2 + 1) * N + idx) * Nfreq + i] = + imag(EH1[i * 6 + k]); + } + } + } + } + + delete[] EH1; + sum_to_master(EH_, EH, 6*2*N*Nfreq); + delete[] EH_; + + /* collapse trailing singleton dimensions */ + while (rank > 0 && dims[rank-1] == 1) + --rank; + /* frequencies are the last dimension */ + if (Nfreq > 1) + dims[++rank] = Nfreq; + + /* output to a file with one dataset per component & real/imag part */ + if (am_master()) { + const int buflen = 1024; + static char filename[buflen]; + snprintf(filename, buflen, "%s%s%s.h5", + prefix ? prefix : "", prefix && prefix[0] ? "-" : "", + fname); + h5file ff(filename, h5file::WRITE, false); + component c[6] = {Ex,Ey,Ez,Hx,Hy,Hz}; + char dataname[128]; + for (int k = 0; k < 6; ++k) + for (int reim = 0; reim < 2; ++reim) { + snprintf(dataname, 128, "%s.%c", + component_name(c[k]), "ri"[reim]); + ff.write(dataname, rank, dims, EH + (k*2 + reim)*N*Nfreq); + } + } + + delete[] EH; +} + +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){ + dft_chunk *F = 0; /* E and H chunks*/ + double eps = 0, mu = 0; + + for (const volume_list *w = where; w; w = w->next) { + direction nd = component_direction(w->c); + if (nd == NO_DIRECTION) nd = normal_direction(w->v); + if (nd == NO_DIRECTION) abort("unknown dft_near2far normal"); + direction fd[2]; + + double weps = get_eps(w->v.center()); + double wmu = get_mu(w->v.center()); + if (w != where && !(approxeq(eps, weps) && approxeq(mu, wmu))) + abort("dft_near2far requires surfaces in a homogeneous medium"); + eps = weps; + mu = wmu; + + /* two transverse directions to normal (in cyclic order to get + correct sign s below) */ + switch (nd) { + case X: fd[0] = Y; fd[1] = Z; break; + case Y: fd[0] = Z; fd[1] = X; break; + case R: fd[0] = P; fd[1] = Z; break; + case P: fd[0] = Z; fd[1] = R; break; + case Z: + if (gv.dim == Dcyl) + fd[0] = R, fd[1] = P; + else + fd[0] = X, fd[1] = Y; + break; + default: abort("invalid normal direction in dft_near2far!"); + } + + 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]); + + /* find equivalent source component c0 and sign s */ + component c0 = direction_component(i == 0 ? Hx : Ex, fd[1-j]); + double s = j == 0 ? 1 : -1; /* sign of n x c */ + if (is_electric(c)) s = -s; + + F = add_dft(c, w->v, freq_min, freq_max, Nfreq, + true, s*w->weight, F, false, 1.0, false, c0); + } + } + } + + return dft_near2far(F, freq_min, freq_max, Nfreq, eps, mu); +} + +} // namespace meep diff --git a/src/stress.cpp b/src/stress.cpp index 4300cd613..8ffd79d01 100644 --- a/src/stress.cpp +++ b/src/stress.cpp @@ -134,10 +134,6 @@ void dft_force::scale_dfts(complex scale) { dft_force fields::add_dft_force(const volume_list *where_, double freq_min, double freq_max, int Nfreq){ dft_chunk *offdiag1 = 0, *offdiag2 = 0, *diag = 0; - direction field_d[3]; - - for (int p = 0; p < 3; ++p) - field_d[p] = gv.yucky_direction(p); volume_list *where = S.reduce(where_); volume_list *where_save = where; diff --git a/tests/Makefile.am b/tests/Makefile.am index 0fb332a36..e03ce8fb2 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -1,8 +1,8 @@ SRC = aniso_disp.cpp bench.cpp bragg_transmission.cpp \ convergence_cyl_waveguide.cpp cylindrical.cpp flux.cpp harmonics.cpp \ -integrate.cpp known_results.cpp one_dimensional.cpp physical.cpp \ -stress_tensor.cpp symmetry.cpp three_d.cpp two_dimensional.cpp \ -2D_convergence.cpp h5test.cpp pml.cpp +integrate.cpp known_results.cpp near2far.cpp one_dimensional.cpp \ +physical.cpp stress_tensor.cpp symmetry.cpp three_d.cpp \ +two_dimensional.cpp 2D_convergence.cpp h5test.cpp pml.cpp EXTRA_DIST = $(SRC) @@ -15,7 +15,7 @@ AM_CPPFLAGS = -I$(top_srcdir)/src .SUFFIXES = .dac .done -check_PROGRAMS = aniso_disp bench bragg_transmission convergence_cyl_waveguide cylindrical flux harmonics integrate known_results one_dimensional physical stress_tensor symmetry three_d two_dimensional 2D_convergence h5test pml +check_PROGRAMS = aniso_disp bench bragg_transmission convergence_cyl_waveguide cylindrical flux harmonics integrate known_results near2far one_dimensional physical stress_tensor symmetry three_d two_dimensional 2D_convergence h5test pml aniso_disp_SOURCES = aniso_disp.cpp aniso_disp_LDADD = $(LIBMEEP) @@ -44,6 +44,9 @@ integrate_LDADD = $(LIBMEEP) known_results_SOURCES = known_results.cpp known_results_LDADD = $(LIBMEEP) +near2far_SOURCES = near2far.cpp +near2far_LDADD = $(LIBMEEP) + one_dimensional_SOURCES = one_dimensional.cpp one_dimensional_LDADD = $(LIBMEEP) @@ -71,7 +74,7 @@ h5test_LDADD = $(LIBMEEP) pml_SOURCES = pml.cpp pml_LDADD = $(LIBMEEP) -TESTS = aniso_disp bench bragg_transmission convergence_cyl_waveguide cylindrical flux harmonics integrate known_results one_dimensional physical stress_tensor symmetry three_d two_dimensional 2D_convergence h5test pml +TESTS = aniso_disp bench bragg_transmission convergence_cyl_waveguide cylindrical flux harmonics integrate known_results near2far one_dimensional physical stress_tensor symmetry three_d two_dimensional 2D_convergence h5test pml LOG_COMPILER = $(RUNCODE) diff --git a/tests/near2far.cpp b/tests/near2far.cpp new file mode 100644 index 000000000..f27bcb9b2 --- /dev/null +++ b/tests/near2far.cpp @@ -0,0 +1,157 @@ +/* Copyright (C) 2005-2014 Massachusetts Institute of Technology +% +% This program is free software; you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation; either version 2, or (at your option) +% any later version. +% +% This program is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with this program; if not, write to the Free Software Foundation, +% Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. +*/ + +/* Check of Green's functions (analytical vs. numerical) + and near-to-far-field transformation. */ + +#include +#include + +#include +using namespace meep; +using namespace std; + +double two(const vec &) { return 2.0; } + +const int EHcomp[10] = {0,1,0,1,2, 3,4,3,4,5}; + +int check_2d_3d(ndim dim, const double xmax, double a, component c0) { + const double dpml = 1; + if (dim != D2 && dim != D3) abort("2d or 3d required"); + grid_volume gv = dim == D2 ? vol2d(xmax + 2*dpml,xmax + 2*dpml,a) : vol3d(xmax + 2*dpml,xmax + 2*dpml,xmax + 2*dpml,a); + gv.center_origin(); + + if (!gv.has_field(c0)) return 1; + master_printf("TESTING %s AT RESOLUTION %g FOR %s SOURCE...\n", + dim == D2 ? "2D" : "3D", a, component_name(c0)); + + structure s(gv, two, pml(dpml)); + fields f(&s); + double w = 0.30; + continuous_src_time src(w); + f.add_point_source(c0, src, zero_vec(dim)); + f.solve_cw(1e-6); + + FOR_E_AND_H(c) if (gv.has_field(c)) { + const int N = 20; + double dx = 0.75 * (xmax/4) / N; + complex F[N], F0[N], EH[6]; + double diff = 0.0, dot0 = 0.0, dot = 0.0; + complex phase = polar(1.0, (4*w*f.dt)*pi); + vec x0 = zero_vec(dim); + for (int i = 0; i < N; ++i) { + double s = xmax/4 + dx*i; + vec x = dim == D2 ? vec(s, 0.5*s) : vec(s, 0.5*s, 0.3*s); + F[i] = f.get_field(c, x) * phase; + (dim == D2 ? green2d : green3d)(EH, x, w, 2.0, 1.0, x0, c0, 1.0); + F0[i] = EH[EHcomp[c]]; + double d = abs(F0[i] - F[i]); + double f = abs(F[i]); + double f0 = abs(F0[i]); + diff += d*d; + dot += f*f; + dot0 += f0*f0; + } + if (dot0 == 0) continue; /* zero field component */ + double relerr = sqrt(diff) / sqrt(dot0); + + master_printf(" GREEN: %s -> %s, resolution %g: relerr = %g\n", + component_name(c0), component_name(c), a, relerr); + + if (relerr > 0.05 * 30/a) { + for (int i = 0; i < N; ++i) + master_printf("%g, %g,%g, %g,%g\n", + xmax/4 + dx*i, + real(F[i]), imag(F[i]), + real(F0[i]), imag(F0[i])); + return 0; + } + } + + const double L = xmax/4; + volume_list vl = + dim == D2 ? + volume_list( volume(vec(-L,+L),vec(+L,+L)), Sy, 1.0, + new volume_list(volume(vec(+L,+L),vec(+L,-L)), Sx, 1.0, + new volume_list(volume(vec(-L,-L),vec(+L,-L)), Sy, -1.0, + new volume_list(volume(vec(-L,-L),vec(-L,+L)), Sx, -1.0)))) : + volume_list( volume(vec(+L,-L,-L),vec(+L,+L,+L)), Sx, +1., + new volume_list(volume(vec(-L,-L,-L),vec(-L,+L,+L)), Sx, -1., + new volume_list(volume(vec(-L,+L,-L),vec(+L,+L,+L)), Sy, +1., + new volume_list(volume(vec(-L,-L,-L),vec(+L,-L,+L)), Sy, -1., + new volume_list(volume(vec(-L,-L,+L),vec(+L,+L,+L)), Sz, +1., + new volume_list(volume(vec(-L,-L,-L),vec(+L,+L,-L)), Sz, -1. + )))))); + + dft_near2far n2f = f.add_dft_near2far(&vl, w, w, 1); + f.update_dfts(); + n2f.scale_dfts(sqrt(2*pi)/f.dt); // cancel time-integration factor + + FOR_E_AND_H(c) if (gv.has_field(c)) { + const int N = 20; + double dx = 0.75 * (xmax/4) / N; + complex F[N], F0[N], EH_[6], EH[6]; + double diff = 0.0, dot = 0.0; + complex phase = polar(1.0, (4*w*f.dt)*pi); + vec x0 = zero_vec(dim); + for (int i = 0; i < N; ++i) { + double s = xmax + dx*i; + vec x = dim == D2 ? vec(s, 0.5*s) : vec(s, 0.5*s, 0.3*s); + n2f.farfield_lowlevel(EH_, x); + sum_to_all(EH_, EH, 6); + F[i] = EH[EHcomp[c]] * phase; + (dim == D2 ? green2d : green3d)(EH, x, w, 2.0, 1.0, x0, c0, 1.0); + F0[i] = EH[EHcomp[c]]; + double d = abs(F0[i] - F[i]); + double f0 = abs(F0[i]); + diff += d*d; + dot += f0*f0; + } + if (dot == 0) continue; /* zero field component */ + double relerr = sqrt(diff) / sqrt(dot); + + master_printf(" NEAR2FAR: %s -> %s, resolution %g: relerr = %g\n", + component_name(c0), component_name(c), a, relerr); + + if (relerr > 0.05 * 30/a) { + for (int i = 0; i < N; ++i) + master_printf("%g, %g,%g, %g,%g\n", + xmax + dx*i, + real(F[i]), imag(F[i]), + real(F0[i]), imag(F0[i])); + return 0; + } + } + + return 1; +} + +int main(int argc, char **argv) { + initialize mpi(argc, argv); + + const double a2d = argc > 1 ? atof(argv[1]) : 20, a3d = argc > 1 ? a2d : 10; + +#if 0 + FOR_E_AND_H(c0) if (!check_2d_3d(D3, 4, a3d, c0)) return 1; +#else + if (!check_2d_3d(D3, 4, a3d, Ez)) return 1; + if (!check_2d_3d(D3, 4, a3d, Hz)) return 1; +#endif + FOR_E_AND_H(c0) if (!check_2d_3d(D2, 8, a2d, c0)) return 1; + + return 0; +}