Near-to-far-field transformation

.gitignore

@@ -73,6 +73,7 @@ tests/h5test
 tests/harmonics
 tests/integrate
 tests/known_results
+tests/near2far
 tests/one_dimensional
 tests/physical
 tests/pml

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<double> (*eps_func)(std::complex<double> omega) = 0,

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,

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

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)		\

src/cw_fields.cpp

@@ -116,6 +116,7 @@ bool fields::solve_cw(double tol, int maxiters, complex<double> 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<double> frequency,
   complex<realnum> *x = reinterpret_cast<complex<realnum>*>(work + nwork);
   complex<realnum> *b = reinterpret_cast<complex<realnum>*>(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

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<double> 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<double> 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<double> weight, dft_chunk *chunk_next,
 			   bool sqrt_dV_and_interp_weights,
 			   complex<double> 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 : 

src/meep.hpp

@@ -827,6 +827,7 @@ class dft_chunk {
 	    std::complex<double> extra_weight_,
 	    component c_,
 	    bool use_centered_grid,
+            ivec shift_, const symmetry &S_, int sn_, int vc,
 	    const void *data_);
   ~dft_chunk();
 
@@ -857,19 +858,22 @@ class dft_chunk {
      it is used in computations involving dft[...], it needs to be public. */
      std::complex<double> extra_weight;
 
-private:
   // parameters passed from field_integrate:
   fields_chunk *fc;
   ivec is, ie;
   vec s0, s1, e0, e1;
   double dV0, dV1;
   bool sqrt_dV_and_interp_weights;
   std::complex<double> scale; // scale factor * phase from shift and symmetry
+  ivec shift;
+  symmetry S; int sn;
 
   // cache of exp(iwt) * scale, of length Nomega
   std::complex<realnum> *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<double> *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<double> *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<double> 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<double> weight = 1.0, dft_chunk *chunk_next = 0,
 		     bool sqrt_dV_and_interp_weights = false,
 		     std::complex<double> 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<double> *EH, const vec &x,
+             double freq, double eps, double mu,
+             const vec &x0, component c0, std::complex<double> f0);
+void green3d(std::complex<double> *EH, const vec &x,
+             double freq, double eps, double mu,
+             const vec &x0, component c0, std::complex<double> f0);
+
 } /* namespace meep */
 
 #endif /* MEEP_H */

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<float> *in, std::complex<double> *out, int size);
 void sum_to_all(const std::complex<double> *in, std::complex<double> *out, int size);
+void sum_to_master(const std::complex<double> *in, std::complex<double> *out, int size);
 long double sum_to_all(long double);
 std::complex<double> sum_to_all(std::complex<double> in);
 std::complex<long double> sum_to_all(std::complex<long double> in);

src/meep/vec.hpp

@@ -918,6 +918,14 @@ class volume_list {
 public:
   volume_list(const volume &v, int c, std::complex<double> 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)

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<float> *in, complex<double> *out, int size) {
   sum_to_all((const float*) in, (double*) out, 2*size);
 }
 
+void sum_to_master(const complex<double> *in, complex<double> *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