Skip to content

Commit

Permalink
updates to ensure correct calculation of eigenmode coefficients with …
Browse files Browse the repository at this point in the history 
…symmetries (NanoComp#417)

* updates to ensure correct calculation of eigenmode coefficients in the presence of symmetries

* store use_symmetry flag in dft_flux structure and abort from get_eigenmode_coefficients() if set incorrectly

* Update Makefile.am
  • Loading branch information
@stevengj
Homer Reid authored and stevengj committed Jul 13, 2018
1 parent cbab237 commit 5c2b12e
Show file tree
Hide file tree
Showing 3 changed files with 37 additions and 27 deletions.
38 changes: 21 additions & 17 deletions src/dft.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -348,26 +348,23 @@ void load_dft_hdf5(dft_chunk *dft_chunks, component c, h5file *file,
load_dft_hdf5(dft_chunks, component_name(c), file, dprefix);
}

dft_flux::dft_flux(const component cE_, const component cH_,
dft_chunk *E_, dft_chunk *H_,
double fmin, double fmax, int Nf,
const volume &where_,
direction normal_direction_) : where(where_)
dft_flux::dft_flux(const component cE_, const component cH_, dft_chunk *E_, dft_chunk *H_,
double fmin, double fmax, int Nf, const volume &where_,
direction normal_direction_, bool use_symmetry_) :
Nfreq(Nf), E(E_), H(H_), cE(cE_), cH(cH_), where(where_),
normal_direction(normal_direction_), use_symmetry(use_symmetry_)
{
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_;
cE = cE_; cH = cH_;
normal_direction = normal_direction_;
}

dft_flux::dft_flux(const dft_flux &f) : where(f.where) {
freq_min = f.freq_min; Nfreq = f.Nfreq; dfreq = f.dfreq;
E = f.E; H = f.H;
cE = f.cE; cH = f.cH;
normal_direction = f.normal_direction;
use_symmetry=f.use_symmetry;
}

double *dft_flux::flux() {
Expand Down Expand Up @@ -417,9 +414,11 @@ void dft_flux::scale_dfts(complex<double> scale) {
}

dft_flux fields::add_dft_flux(const volume_list *where_,
double freq_min, double freq_max, int Nfreq) {
double freq_min, double freq_max, int Nfreq,
bool use_symmetry)
{
if (!where_) // handle empty list of volumes
return dft_flux(Ex, Hy, NULL, NULL, freq_min, freq_max, Nfreq, v, NO_DIRECTION);
return dft_flux(Ex, Hy, NULL, NULL, freq_min, freq_max, Nfreq, v, NO_DIRECTION, use_symmetry);

dft_chunk *E = 0, *H = 0;
component cE[2] = {Ex,Ey}, cH[2] = {Hy,Hx};
Expand All @@ -430,7 +429,7 @@ dft_flux fields::add_dft_flux(const volume_list *where_,
// to store the first volume in the list.
volume firstvol(where_->v);

volume_list *where = S.reduce(where_);
volume_list *where = use_symmetry ? S.reduce(where_) : new volume_list(where_);
volume_list *where_save = where;
while (where) {
derived_component c = derived_component(where->c);
Expand Down Expand Up @@ -465,7 +464,7 @@ dft_flux fields::add_dft_flux(const volume_list *where_,
// if the volume list has only one entry, store its component's direction.
// if the volume list has > 1 entry, store NO_DIRECTION.
direction flux_dir = (where_->next ? NO_DIRECTION : component_direction(where_->c));
return dft_flux(cE[0], cH[0], E, H, freq_min, freq_max, Nfreq, firstvol, flux_dir);
return dft_flux(cE[0], cH[0], E, H, freq_min, freq_max, Nfreq, firstvol, flux_dir, use_symmetry);
}


Expand All @@ -486,15 +485,20 @@ direction fields::normal_direction(const volume &where) const {
}

dft_flux fields::add_dft_flux(direction d, const volume &where,
double freq_min, double freq_max, int Nfreq) {
double freq_min, double freq_max, int Nfreq,
bool use_symmetry) {
if (d == NO_DIRECTION)
d = normal_direction(where);
volume_list vl(where, direction_component(Sx, d));
dft_flux flux=add_dft_flux(&vl, freq_min, freq_max, Nfreq);
dft_flux flux=add_dft_flux(&vl, freq_min, freq_max, Nfreq, use_symmetry);
flux.normal_direction=d;
return flux;
}

dft_flux fields::add_mode_monitor(direction d, const volume &where,
double freq_min, double freq_max, int Nfreq)
{ return add_dft_flux(d,where,freq_min,freq_max,Nfreq, /*use_symmetry=*/false); }

dft_flux fields::add_dft_flux_box(const volume &where,
double freq_min, double freq_max, int Nfreq){
volume_list *faces = 0;
Expand Down Expand Up @@ -646,9 +650,9 @@ cdouble dft_chunk::process_dft_component(int rank, direction *ds,

cdouble mode1val=0.0, mode2val=0.0;
if (mode1_data)
mode1val=eigenmode_amplitude(mode1_data,loc,c_conjugate);
mode1val=eigenmode_amplitude(mode1_data,loc,S.transform(c_conjugate,sn));
if (mode2_data)
mode2val=eigenmode_amplitude(mode2_data,loc,c);
mode2val=eigenmode_amplitude(mode2_data,loc,S.transform(c,sn));

if (file)
{ int idx2 = ((((file_offset[0] + file_offset[1] + file_offset[2])
Expand Down
14 changes: 10 additions & 4 deletions src/meep.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -935,7 +935,8 @@ class dft_flux {
dft_flux(const component cE_, const component cH_,
dft_chunk *E_, dft_chunk *H_,
double fmin, double fmax, int Nf,
const volume &where_, direction normal_direction_);
const volume &where_, direction normal_direction_,
bool use_symmetry_);
dft_flux(const dft_flux &f);

double *flux();
Expand All @@ -960,6 +961,7 @@ class dft_flux {
component cE, cH;
volume where;
direction normal_direction;
bool use_symmetry;
};

// stress.cpp (normally created with fields::add_dft_force)
Expand Down Expand Up @@ -1539,14 +1541,18 @@ class fields {
double freq_min, double freq_max, int Nfreq,
bool include_dV = true);
void update_dfts();
dft_flux add_dft_flux(const volume_list *where,
double freq_min, double freq_max, int Nfreq, bool use_symmetry=true);
dft_flux add_dft_flux(direction d, const volume &where,
double freq_min, double freq_max, int Nfreq);
double freq_min, double freq_max, int Nfreq, bool use_symmetry=true);
dft_flux add_dft_flux_box(const volume &where,
double freq_min, double freq_max, int Nfreq);
dft_flux add_dft_flux_plane(const volume &where,
double freq_min, double freq_max, int Nfreq);
dft_flux add_dft_flux(const volume_list *where,
double freq_min, double freq_max, int Nfreq);

// a "mode monitor" is just a dft_flux with symmetry reduction turned off.
dft_flux add_mode_monitor(direction d, const volume &where,
double freq_min, double freq_max, int Nfreq);

dft_fields add_dft_fields(component *components, int num_components,
const volume where,
Expand Down
12 changes: 6 additions & 6 deletions src/mpb.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -230,10 +230,10 @@ void *fields::get_eigenmode(double omega_src,
// in any direction, set the corresponding component of the eigenmode initial-guess
// k-vector to be the (real part of the) bloch vector in that direction.
vec kpoint(_kpoint);
LOOP_OVER_DIRECTIONS(v.dim, d)
if (float(eig_vol.in_direction(d) == float(v.in_direction(d))))
if (boundaries[High][d]==Periodic && boundaries[Low][d]==Periodic)
kpoint.set_direction(d, real(k[d]));
LOOP_OVER_DIRECTIONS(v.dim, dd)
if (float(eig_vol.in_direction(dd) == float(v.in_direction(dd))))
if (boundaries[High][dd]==Periodic && boundaries[Low][dd]==Periodic)
kpoint.set_direction(dd, real(k[dd]));

//bool verbose=true;
if (resolution <= 0.0) resolution = 2 * gv.a; // default to twice resolution
Expand Down Expand Up @@ -634,8 +634,8 @@ void fields::get_eigenmode_coefficients(dft_flux flux,
if (d==NO_DIRECTION)
abort("cannot determine normal direction in get_eigenmode_coefficients");

if (S.multiplicity() > 1)
abort("symmetries are not yet supported in get_eigenmode_coefficients");
if (flux.use_symmetry && S.multiplicity() > 1)
abort("flux regions for eigenmode projection must be created by add_mode_monitor()");

vec kpoint(0.0,0.0,0.0); // default guess

Expand Down

0 comments on commit 5c2b12e

Please sign in to comment.