From b40c8087d4f04db8289ed57938db8e8cc793a74e Mon Sep 17 00:00:00 2001 From: Ardavan Oskooi Date: Wed, 27 Feb 2019 08:35:26 -0800 Subject: [PATCH] dft_energy routines and Scheme SWIG wrappers (#744) --- scheme/meep-ctl-swig.hpp | 3 + scheme/meep.cpp | 24 +++++++ scheme/meep.scm.in | 68 ++++++++++++++++++ src/dft.cpp | 144 +++++++++++++++++++++++++++++++++++++++ src/meep.hpp | 37 ++++++++++ 5 files changed, 276 insertions(+) diff --git a/scheme/meep-ctl-swig.hpp b/scheme/meep-ctl-swig.hpp index fcee63482..471068ee6 100644 --- a/scheme/meep-ctl-swig.hpp +++ b/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); diff --git a/scheme/meep.cpp b/scheme/meep.cpp index e4359a5cd..61f5e5980 100644 --- a/scheme/meep.cpp +++ b/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; diff --git a/scheme/meep.scm.in b/scheme/meep.scm.in index 90eb61fc9..fee50d199 100644 --- a/scheme/meep.scm.in +++ b/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 diff --git a/src/dft.cpp b/src/dft.cpp index fe1527f90..184979e11 100644 --- a/src/dft.cpp +++ b/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 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) { diff --git a/src/meep.hpp b/src/meep.hpp index 38048f547..bdc2cb6ee 100644 --- a/src/meep.hpp +++ b/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 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 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);