dft_energy routines and Scheme SWIG wrappers (#744)

scheme/meep-ctl-swig.hpp

@@ -38,6 +38,9 @@ kpoint_list do_get_eigenmode_coefficients(meep::fields *f, meep::dft_flux flux,
                                           void *user_kpoint_data);
 
 ctlio::number_list dft_flux_flux(meep::dft_flux *f);
+ctlio::number_list dft_energy_electric(meep::dft_energy *f);
+ctlio::number_list dft_energy_magnetic(meep::dft_energy *f);
+ctlio::number_list dft_energy_total(meep::dft_energy *f);
 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);

scheme/meep.cpp

@@ -95,6 +95,30 @@ ctlio::number_list dft_flux_flux(dft_flux *f) {
   return res;
 }
 
+ctlio::number_list dft_energy_electric(dft_energy *f)
+{
+  ctlio::number_list res;
+  res.num_items = f->Nfreq;
+  res.items = f->electric();
+  return res;
+}
+
+ctlio::number_list dft_energy_magnetic(dft_energy *f)
+{
+  ctlio::number_list res;
+  res.num_items = f->Nfreq;
+  res.items = f->magnetic();
+  return res;
+}
+
+ctlio::number_list dft_energy_total(dft_energy *f)
+{
+  ctlio::number_list res;
+  res.num_items = f->Nfreq;
+  res.items = f->total();
+  return res;
+}
+
 ctlio::number_list dft_force_force(dft_force *f) {
   ctlio::number_list res;
   res.num_items = f->Nfreq;

scheme/meep.scm.in

@@ -604,6 +604,74 @@
   (if (null? fields) (init-fields))
   (apply fields-add-mode-monitor (append (list fields fcen df nfreq) fluxes)))
 
+; ****************************************************************
+; Energy spectra
+
+(define-class energy-region no-parent
+  (define-property center no-default 'vector3)
+  (define-property size (vector3 0 0 0) 'vector3)
+  (define-property direction 0 'integer)
+  (define-property weight 1.0 'cnumber))
+
+(define (fields-add-energy fields fcen df nfreq . energys)
+  (fields-add-fluxish-stuff meep-fields-add-dft-energy
+			    fields fcen df nfreq energys))
+
+(define (add-energy fcen df nfreq . energys)
+  (if (null? fields) (init-fields))
+  (apply fields-add-energy (append (list fields fcen df nfreq) energys)))
+
+(define (scale-energy-fields s f)
+  (meep-dft-energy-scale-dfts f s))
+
+(define (get-energy-freqs f)
+  (arith-sequence
+   (meep-dft-energy-freq-min-get f)
+   (meep-dft-energy-dfreq-get f)
+   (meep-dft-energy-Nfreq-get f)))
+
+(define (get-electric-energy f)
+  (dft-energy-electric f))
+
+(define (get-magnetic-energy f)
+  (dft-energy-magnetic f))
+
+(define (get-total-energy f)
+  (dft-energy-total f))
+
+(define (display-electric-energy . energys)
+  (if (not (null? energys))
+      (apply display-csv
+	     (append (list "electric-energy"
+			   (get-energy-freqs (car energys)))
+		     (map get-electric-energy energys)))))
+
+(define (display-magnetic-energy . energys)
+  (if (not (null? energys))
+      (apply display-csv
+	     (append (list "magnetic-energy"
+			   (get-energy-freqs (car energys)))
+		     (map get-magnetic-energy energys)))))
+
+(define (display-total-energy . energys)
+  (if (not (null? energys))
+      (apply display-csv
+	     (append (list "total-energy"
+			   (get-energy-freqs (car energys)))
+		     (map get-total-energy energys)))))
+
+(define (load-energy fname energy)
+  (if (null? fields) (init-fields))
+  (meep-dft-energy-load-hdf5 energy fields fname "" (get-filename-prefix)))
+
+(define (save-energy fname energy)
+  (if (null? fields) (init-fields))
+  (meep-dft-energy-save-hdf5 energy fields fname "" (get-filename-prefix)))
+
+(define (load-minus-energy fname energy)
+  (load-energy fname energy)
+  (meep-dft-energy-scale-dfts energy -1.0))
+
 ; ****************************************************************
 ; Force spectra (from stress tensor) - very similar interface to flux spectra
 

src/dft.cpp

@@ -444,6 +444,150 @@ dft_flux fields::add_dft_flux(const volume_list *where_, double freq_min, double
   return dft_flux(cE[0], cH[0], E, H, freq_min, freq_max, Nfreq, firstvol, flux_dir, use_symmetry);
 }
 
+dft_energy::dft_energy(dft_chunk *E_, dft_chunk *H_,
+                       dft_chunk *D_, dft_chunk *B_,
+                       double fmin, double fmax, int Nf)
+{
+  if (Nf <= 1) fmin = fmax = (fmin + fmax) * 0.5;
+  freq_min = fmin;
+  Nfreq = Nf;
+  dfreq = Nf <= 1 ? 0.0 : (fmax - fmin) / (Nf - 1);
+  E = E_; H = H_; D = D_; B = B_;
+}
+
+dft_energy::dft_energy(const dft_energy &f) {
+  freq_min = f.freq_min; Nfreq = f.Nfreq; dfreq = f.dfreq;
+  E = f.E; H = f.H; D = f.D; B = f.B;
+}
+
+double *dft_energy::electric() {
+  double *F = new double[Nfreq];
+  for (int i = 0; i < Nfreq; ++i) F[i] = 0;
+  for (dft_chunk *curE = E, *curD = D; curE && curD;
+       curE = curE->next_in_dft, curD = curD->next_in_dft)
+    for (size_t k = 0; k < curE->N; ++k)
+      for (int i = 0; i < Nfreq; ++i)
+	F[i] += 0.5*real(conj(curE->dft[k*Nfreq + i])
+			 * curD->dft[k*Nfreq + i]);
+  double *Fsum = new double[Nfreq];
+  sum_to_all(F, Fsum, Nfreq);
+  delete[] F;
+  return Fsum;
+}
+
+double *dft_energy::magnetic() {
+  double *F = new double[Nfreq];
+  for (int i = 0; i < Nfreq; ++i) F[i] = 0;
+  for (dft_chunk *curH = H, *curB = B; curH && curB;
+       curH = curH->next_in_dft, curB = curB->next_in_dft)
+    for (size_t k = 0; k < curH->N; ++k)
+      for (int i = 0; i < Nfreq; ++i)
+	F[i] += 0.5*real(conj(curH->dft[k*Nfreq + i])
+			 * curB->dft[k*Nfreq + i]);
+  double *Fsum = new double[Nfreq];
+  sum_to_all(F, Fsum, Nfreq);
+  delete[] F;
+  return Fsum;
+}
+
+double *dft_energy::total() {
+  double *Fe = electric();
+  double *Fm = magnetic();
+  double *F = new double[Nfreq];
+  for (int i = 0; i < Nfreq; ++i) F[i] = Fe[i]+Fm[i];
+  delete[] Fe;
+  delete[] Fm;
+  return F;
+}
+
+dft_energy fields::add_dft_energy(const volume_list *where_,
+                                  double freq_min, double freq_max, int Nfreq) {
+  dft_chunk *E = 0, *D = 0, *H = 0, *B = 0;
+  volume_list *where = S.reduce(where_);
+  volume_list *where_save = where;
+  while (where) {
+    LOOP_OVER_FIELD_DIRECTIONS(gv.dim, d) {
+      E = add_dft(direction_component(Ex, d), where->v, freq_min, freq_max, Nfreq,
+		  true, 1.0, E);
+      D = add_dft(direction_component(Dx, d), where->v, freq_min, freq_max, Nfreq,
+		  true, 1.0, D);
+      H = add_dft(direction_component(Hx, d), where->v, freq_min, freq_max, Nfreq,
+		  true, 1.0, H);
+      B = add_dft(direction_component(Bx, d), where->v, freq_min, freq_max, Nfreq,
+		  true, 1.0, B);
+    }
+    where = where->next;
+  }
+  delete where_save;
+
+  return dft_energy(E, H, D, B, freq_min, freq_max, Nfreq);
+}
+
+void dft_energy::save_hdf5(h5file *file, const char *dprefix) {
+  save_dft_hdf5(E, "E", file, dprefix);
+  file->prevent_deadlock(); // hackery
+  save_dft_hdf5(D, "D", file, dprefix);
+  file->prevent_deadlock(); // hackery
+  save_dft_hdf5(H, "H", file, dprefix);
+  file->prevent_deadlock(); // hackery
+  save_dft_hdf5(B, "B", file, dprefix);
+}
+
+void dft_energy::load_hdf5(h5file *file, const char *dprefix) {
+  load_dft_hdf5(E, "E", file, dprefix);
+  file->prevent_deadlock(); // hackery
+  load_dft_hdf5(D, "D", file, dprefix);
+  file->prevent_deadlock(); // hackery
+  load_dft_hdf5(H, "H", file, dprefix);
+  file->prevent_deadlock(); // hackery
+  load_dft_hdf5(B, "B", file, dprefix);
+}
+
+void dft_energy::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_energy::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_energy::scale_dfts(complex<double> scale) {
+  if (E) E->scale_dft(scale);
+  if (D) D->scale_dft(scale);
+  if (H) H->scale_dft(scale);
+  if (B) B->scale_dft(scale);
+}
+
+void dft_energy::remove()
+{
+  while (E) {
+    dft_chunk *nxt = E->next_in_dft;
+    delete E;
+    E = nxt;
+  }
+  while (D) {
+    dft_chunk *nxt = D->next_in_dft;
+    delete D;
+    D = nxt;
+  }
+  while (H) {
+    dft_chunk *nxt = H->next_in_dft;
+    delete H;
+    H = nxt;
+  }
+  while (B) {
+    dft_chunk *nxt = B->next_in_dft;
+    delete B;
+    B = nxt;
+  }
+}
+
 direction fields::normal_direction(const volume &where) const {
   direction d = where.normal_direction();
   if (d == NO_DIRECTION) {

src/meep.hpp

@@ -1049,6 +1049,41 @@ class dft_flux {
   bool use_symmetry;
 };
 
+// dft.cpp (normally created with fields::add_dft_energy)
+class dft_energy {
+public:
+  dft_energy(dft_chunk *E_, dft_chunk *H_, dft_chunk *D_, dft_chunk *B_,
+             double fmin, double fmax, int Nf);
+  dft_energy(const dft_energy &f);
+
+  double *electric();
+  double *magnetic();
+  double *total();
+
+  void save_hdf5(h5file *file, const char *dprefix = 0);
+  void load_hdf5(h5file *file, const char *dprefix = 0);
+
+  void operator-=(const dft_energy &fl) {
+    if (E && fl.E) *E -= *fl.E;
+    if (H && fl.H) *H -= *fl.H;
+    if (D && fl.D) *D -= *fl.D;
+    if (B && fl.B) *B -= *fl.B;
+  }
+
+  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 *E, *H, *D, *B;
+};
+
 // stress.cpp (normally created with fields::add_dft_force)
 class dft_force {
 public:
@@ -1687,6 +1722,8 @@ class fields {
   void get_mode_mode_overlap(void *mode1_data, void *mode2_data, dft_flux flux,
                              std::complex<double> overlaps[2]);
 
+  dft_energy add_dft_energy(const volume_list *where, double freq_min, double freq_max, int Nfreq);
+
   // stress.cpp
   dft_force add_dft_force(const volume_list *where, double freq_min, double freq_max, int Nfreq);