diff --git a/scheme/examples/group-velocity.ctl b/scheme/examples/group-velocity.ctl new file mode 100644 index 000000000..050777f8c --- /dev/null +++ b/scheme/examples/group-velocity.ctl @@ -0,0 +1,37 @@ +;; compute group velocity of a waveguide mode using two different methods +;; (1) ratio of Poynting flux to energy density +;; (2) via MPB from get-eigenmode-coefficients + +(set-param! resolution 20) + +(set! geometry-lattice (make lattice (size 10 5 no-size))) + +(set! geometry (list (make block + (center 0 0 0) + (size infinity 1 infinity) + (material (make medium (epsilon 12)))))) + +(set! pml-layers (list (make pml (thickness 1)))) + +(define-param fsrc 0.15) + +(set! sources (list (make eigenmode-source + (src (make gaussian-src (frequency fsrc) (fwidth (* 0.2 fsrc)))) + (center -3 0 0) + (size 0 5 0) + (eig-band 1) + (eig-parity (+ ODD-Z EVEN-Y)) + (eig-match-freq? true)))) + +(define flux (add-flux fsrc 0 1 (make flux-region (center 3 0 0) (size 0 5 0)))) +(define energy (add-energy fsrc 0 1 (make energy-region (center 3 0 0) (size 0 5 0)))) +(run-sources+ 100) + +(define res (get-eigenmode-coefficients flux (list 1) #:eig-parity (+ ODD-Z EVEN-Y))) +(define mode-vg (array-ref (list-ref res 1) 0 0)) + +(define poynting-flux (list-ref (get-fluxes flux) 0)) +(define e-energy (list-ref (get-electric-energy energy) 0)) +(define ratio-vg (/ (* 0.5 poynting-flux) e-energy)) + +(print "group-velocity:, " ratio-vg ", " mode-vg "\n") diff --git a/src/dft.cpp b/src/dft.cpp index 184979e11..b44a1ee0c 100644 --- a/src/dft.cpp +++ b/src/dft.cpp @@ -446,29 +446,32 @@ dft_flux fields::add_dft_flux(const volume_list *where_, double freq_min, double dft_energy::dft_energy(dft_chunk *E_, dft_chunk *H_, dft_chunk *D_, dft_chunk *B_, - double fmin, double fmax, int Nf) -{ + double fmin, double fmax, int Nf, const volume &where_) + : Nfreq(Nf), E(E_), H(H_), D(D_), B(B_), where(where_) { 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; +dft_energy::dft_energy(const dft_energy &f) : where(f.where) { + 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 (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]); + 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; @@ -477,13 +480,13 @@ double *dft_energy::electric() { double *dft_energy::magnetic() { double *F = new double[Nfreq]; - for (int i = 0; i < Nfreq; ++i) F[i] = 0; + 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]); + 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; @@ -494,7 +497,8 @@ 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]; + for (int i = 0; i < Nfreq; ++i) + F[i] = Fe[i]+Fm[i]; delete[] Fe; delete[] Fm; return F; @@ -502,25 +506,26 @@ double *dft_energy::total() { dft_energy fields::add_dft_energy(const volume_list *where_, double freq_min, double freq_max, int Nfreq) { + + if (!where_) // handle empty list of volumes + return dft_energy(NULL, NULL, NULL, NULL, freq_min, freq_max, Nfreq, v); + dft_chunk *E = 0, *D = 0, *H = 0, *B = 0; - volume_list *where = S.reduce(where_); + volume firstvol(where_->v); + volume_list *where = new volume_list(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); + 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, false, 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, false, 1.0, B); } where = where->next; } delete where_save; - return dft_energy(E, H, D, B, freq_min, freq_max, Nfreq); + return dft_energy(E, H, D, B, freq_min, freq_max, Nfreq, firstvol); } void dft_energy::save_hdf5(h5file *file, const char *dprefix) { diff --git a/src/meep.hpp b/src/meep.hpp index 8a3e7b9fd..6a1072b11 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -1054,7 +1054,7 @@ class dft_flux { class dft_energy { public: dft_energy(dft_chunk *E_, dft_chunk *H_, dft_chunk *D_, dft_chunk *B_, - double fmin, double fmax, int Nf); + double fmin, double fmax, int Nf, const volume &where_); dft_energy(const dft_energy &f); double *electric(); @@ -1083,6 +1083,7 @@ class dft_energy { double freq_min, dfreq; int Nfreq; dft_chunk *E, *H, *D, *B; + volume where; }; // stress.cpp (normally created with fields::add_dft_force)