Skip to content

Commit

Permalink
Add missing MPB functionality (NanoComp#223)
Browse files Browse the repository at this point in the history 
* Add first_brillouin_zone and test

* Add missing functionality

* Add missing functionality
  • Loading branch information
@ChristopherHogan @stevengj
ChristopherHogan authored and stevengj committed Mar 9, 2018
1 parent 8fc21ff commit 4fdadb4
Show file tree
Hide file tree
Showing 9 changed files with 708 additions and 20 deletions.
361 changes: 351 additions & 10 deletions libpympb/pympb.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,18 @@
int xyz_index = ((i2_ * n1 + i1) * n3 + i3);
# endif /* HAVE_MPI */

#ifdef CASSIGN_MULT
#undef CASSIGN_MULT
#endif

/* a = b * c */
#define CASSIGN_MULT(a, b, c) { \
mpb_real bbbb_re = (b).re, bbbb_im = (b).im; \
mpb_real cccc_re = (c).re, cccc_im = (c).im; \
CASSIGN_SCALAR(a, bbbb_re * cccc_re - bbbb_im * cccc_im, \
bbbb_re * cccc_im + bbbb_im * cccc_re); \
}

// TODO: Support MPI
#define mpi_allreduce(sb, rb, n, ctype, t, op, comm) { \
CHECK((sb) != (rb), "MPI_Allreduce doesn't work for sendbuf == recvbuf");\
Expand Down Expand Up @@ -87,18 +99,20 @@ static int mean_epsilon_func(symmetric_matrix* meps, symmetric_matrix *meps_inv,
/* a couple of utilities to convert libctl data types to the data
types of the eigensolver & maxwell routines: */

void vector3_to_arr(mpb_real arr[3], vector3 v)
{
arr[0] = v.x;
arr[1] = v.y;
arr[2] = v.z;
void vector3_to_arr(mpb_real arr[3], vector3 v) {
arr[0] = v.x;
arr[1] = v.y;
arr[2] = v.z;
}

void matrix3x3_to_arr(mpb_real arr[3][3], matrix3x3 m)
{
vector3_to_arr(arr[0], m.c0);
vector3_to_arr(arr[1], m.c1);
vector3_to_arr(arr[2], m.c2);
void matrix3x3_to_arr(mpb_real arr[3][3], matrix3x3 m) {
vector3_to_arr(arr[0], m.c0);
vector3_to_arr(arr[1], m.c1);
vector3_to_arr(arr[2], m.c2);
}

cnumber cscalar2cnumber(scalar_complex cs) {
return make_cnumber(CSCALAR_RE(cs), CSCALAR_IM(cs));
}

// Return a string describing the current parity, used for frequency and filename
Expand Down Expand Up @@ -1202,6 +1216,69 @@ void mode_solver::get_epsilon() {
eps_high == eps_low ? 100.0 : 100.0 * (eps_mean-eps_low) / (eps_high-eps_low));
}

/* get the mu function, and compute some statistics */
void mode_solver::get_mu() {
mpb_real eps_mean = 0;
mpb_real mu_inv_mean = 0;
mpb_real eps_high = -1e20;
mpb_real eps_low = 1e20;
int fill_count = 0;

if (!mdata) {
meep::master_fprintf(stderr, "mode_solver.init must be called before get-mu!\n");
return;
}

curfield = (scalar_complex *) mdata->fft_data;
mpb_real *mu = (mpb_real *) curfield;
curfield_band = 0;
curfield_type = mu_CURFIELD_TYPE;

/* get mu. Recall that we actually have an inverse
dielectric tensor at each point; define an average index by
the inverse of the average eigenvalue of the 1/eps tensor.
i.e. 3/(trace 1/eps). */

int N = mdata->fft_output_size;

for (int i = 0; i < N; ++i) {
if (mdata->mu_inv == NULL) {
mu[i] = 1.0;
}
else {
mu[i] = mean_medium_from_matrix(mdata->mu_inv + i);
}

if (mu[i] < eps_low) {
eps_low = mu[i];
}
if (mu[i] > eps_high) {
eps_high = mu[i];
}

eps_mean += mu[i];
mu_inv_mean += 1/mu[i];
if (mu[i] > 1.0001) {
++fill_count;
}
}

mpi_allreduce_1(&eps_mean, mpb_real, SCALAR_MPI_TYPE, MPI_SUM, mpb_comm);
mpi_allreduce_1(&mu_inv_mean, mpb_real, SCALAR_MPI_TYPE, MPI_SUM, mpb_comm);
mpi_allreduce_1(&eps_low, mpb_real, SCALAR_MPI_TYPE, MPI_MIN, mpb_comm);
mpi_allreduce_1(&eps_high, mpb_real, SCALAR_MPI_TYPE, MPI_MAX, mpb_comm);
mpi_allreduce_1(&fill_count, int, MPI_INT, MPI_SUM, mpb_comm);

N = mdata->nx * mdata->ny * mdata->nz;
eps_mean /= N;
mu_inv_mean = N/mu_inv_mean;

meep::master_printf("mu: %g-%g, mean %g, harm. mean %g, %g%% > 1, %g%% \"fill\"\n",
eps_low, eps_high, eps_mean, mu_inv_mean, (100.0 * fill_count) / N,
eps_high == eps_low ? 100.0 : 100.0 *
(eps_mean-eps_low) / (eps_high-eps_low));
}

void mode_solver::curfield_reset() {
curfield = NULL;
curfield_type = '-';
Expand Down Expand Up @@ -1595,6 +1672,168 @@ std::vector<mpb_real> mode_solver::get_output_k() {
return output_k;
}

mpb_real mode_solver::get_val(int ix, int iy, int iz, int nx, int ny, int nz,
int last_dim_size, mpb_real *data, int stride, int conjugate) {
// #ifdef HAVE_MPI
// CHECK(0, "get-*-point not yet implemented for MPI!");
// #else
(void)nx;
(void)last_dim_size;
(void)conjugate;
return data[(((ix * ny) + iy) * nz + iz) * stride];
// #endif
}

mpb_real mode_solver::interp_val(vector3 p, int nx, int ny, int nz, int last_dim_size,
mpb_real *data, int stride, int conjugate) {
double ipart;
mpb_real rx, ry, rz, dx, dy, dz;
int x, y, z, x2, y2, z2;

mpb_real latx = geometry_lattice.size.x == 0 ? 1e-20 : geometry_lattice.size.x;
mpb_real laty = geometry_lattice.size.y == 0 ? 1e-20 : geometry_lattice.size.y;
mpb_real latz = geometry_lattice.size.z == 0 ? 1e-20 : geometry_lattice.size.z;

rx = modf(p.x/latx + 0.5, &ipart); if (rx < 0) rx += 1;
ry = modf(p.y/laty + 0.5, &ipart); if (ry < 0) ry += 1;
rz = modf(p.z/latz + 0.5, &ipart); if (rz < 0) rz += 1;

/* get the point corresponding to r in the grid: */
x = rx * nx;
y = ry * ny;
z = rz * nz;

/* get the difference between (x,y,z) and the actual point */
dx = rx * nx - x;
dy = ry * ny - y;
dz = rz * nz - z;

/* get the other closest point in the grid, with periodic boundaries: */
x2 = (nx + (dx >= 0.0 ? x + 1 : x - 1)) % nx;
y2 = (ny + (dy >= 0.0 ? y + 1 : y - 1)) % ny;
z2 = (nz + (dz >= 0.0 ? z + 1 : z - 1)) % nz;

/* take abs(d{xyz}) to get weights for {xyz} and {xyz}2: */
dx = fabs(dx);
dy = fabs(dy);
dz = fabs(dz);

#define D(x,y,z) (get_val(x,y,z,nx,ny,nz,last_dim_size, data,stride,conjugate))

return(((D(x,y,z)*(1.0-dx) + D(x2,y,z)*dx) * (1.0-dy) +
(D(x,y2,z)*(1.0-dx) + D(x2,y2,z)*dx) * dy) * (1.0-dz) +
((D(x,y,z2)*(1.0-dx) + D(x2,y,z2)*dx) * (1.0-dy) +
(D(x,y2,z2)*(1.0-dx) + D(x2,y2,z2)*dx) * dy) * dz);

#undef D
}

scalar_complex mode_solver::interp_cval(vector3 p, int nx, int ny, int nz,
int last_dim_size, mpb_real *data, int stride) {
scalar_complex cval;
cval.re = interp_val(p, nx,ny,nz,last_dim_size, data, stride, 0);
cval.im = interp_val(p, nx,ny,nz,last_dim_size,data + 1, stride, 1);
return cval;
}

#define f_interp_val(p,f,data,stride,conj) interp_val(p,f->nx,f->ny,f->nz,f->last_dim_size,data,stride,conj)
#define f_interp_cval(p,f,data,stride) interp_cval(p,f->nx,f->ny,f->nz,f->last_dim_size,data,stride)

symmetric_matrix mode_solver::interp_eps_inv(vector3 p) {
int stride = sizeof(symmetric_matrix) / sizeof(mpb_real);
symmetric_matrix eps_inv;

eps_inv.m00 = f_interp_val(p, mdata, &mdata->eps_inv->m00, stride, 0);
eps_inv.m11 = f_interp_val(p, mdata, &mdata->eps_inv->m11, stride, 0);
eps_inv.m22 = f_interp_val(p, mdata, &mdata->eps_inv->m22, stride, 0);
#ifdef WITH_HERMITIAN_EPSILON
eps_inv.m01 = f_interp_cval(p, mdata, &mdata->eps_inv->m01.re, stride);
eps_inv.m02 = f_interp_cval(p, mdata, &mdata->eps_inv->m02.re, stride);
eps_inv.m12 = f_interp_cval(p, mdata, &mdata->eps_inv->m12.re, stride);
#else
eps_inv.m01 = f_interp_val(p, mdata, &mdata->eps_inv->m01, stride, 0);
eps_inv.m02 = f_interp_val(p, mdata, &mdata->eps_inv->m02, stride, 0);
eps_inv.m12 = f_interp_val(p, mdata, &mdata->eps_inv->m12, stride, 0);
#endif
return eps_inv;
}

mpb_real mode_solver::get_epsilon_point(vector3 p) {
symmetric_matrix eps_inv;
eps_inv = interp_eps_inv(p);
return mean_medium_from_matrix(&eps_inv);
}

cmatrix3x3 mode_solver::get_epsilon_inverse_tensor_point(vector3 p) {
symmetric_matrix eps_inv;
eps_inv = interp_eps_inv(p);

#ifdef WITH_HERMITIAN_EPSILON
return make_hermitian_cmatrix3x3(eps_inv.m00,eps_inv.m11,eps_inv.m22,
cscalar2cnumber(eps_inv.m01),
cscalar2cnumber(eps_inv.m02),
cscalar2cnumber(eps_inv.m12));
#else
return make_hermitian_cmatrix3x3(eps_inv.m00,eps_inv.m11,eps_inv.m22,
make_cnumber(eps_inv.m01,0),
make_cnumber(eps_inv.m02,0),
make_cnumber(eps_inv.m12,0));
#endif
}

mpb_real mode_solver::get_energy_point(vector3 p) {
CHECK(curfield && strchr("DHBR", curfield_type),
"compute-field-energy must be called before get-energy-point");
return f_interp_val(p, mdata, (mpb_real *) curfield, 1, 0);
}

cvector3 mode_solver::get_bloch_field_point(vector3 p) {
scalar_complex field[3];
cvector3 F;

CHECK(curfield && strchr("dhbecv", curfield_type),
"field must be must be loaded before get-*field*-point");
field[0] = f_interp_cval(p, mdata, &curfield[0].re, 6);
field[1] = f_interp_cval(p, mdata, &curfield[1].re, 6);
field[2] = f_interp_cval(p, mdata, &curfield[2].re, 6);
F.x = cscalar2cnumber(field[0]);
F.y = cscalar2cnumber(field[1]);
F.z = cscalar2cnumber(field[2]);
return F;
}

cvector3 mode_solver::get_field_point(vector3 p) {
scalar_complex field[3], phase;
cvector3 F;

CHECK(curfield && strchr("dhbecv", curfield_type),
"field must be must be loaded before get-*field*-point");
field[0] = f_interp_cval(p, mdata, &curfield[0].re, 6);
field[1] = f_interp_cval(p, mdata, &curfield[1].re, 6);
field[2] = f_interp_cval(p, mdata, &curfield[2].re, 6);

if (curfield_type != 'v') {
mpb_real latx = geometry_lattice.size.x == 0 ? 1e-20 : geometry_lattice.size.x;
mpb_real laty = geometry_lattice.size.y == 0 ? 1e-20 : geometry_lattice.size.y;
mpb_real latz = geometry_lattice.size.z == 0 ? 1e-20 : geometry_lattice.size.z;

double phase_phi = TWOPI * (cur_kvector.x * (p.x / latx) +
cur_kvector.y * (p.y / laty) +
cur_kvector.z * (p.z / latz));

CASSIGN_SCALAR(phase, cos(phase_phi), sin(phase_phi));
CASSIGN_MULT(field[0], field[0], phase);
CASSIGN_MULT(field[1], field[1], phase);
CASSIGN_MULT(field[2], field[2], phase);
}

F.x = cscalar2cnumber(field[0]);
F.y = cscalar2cnumber(field[1]);
F.z = cscalar2cnumber(field[2]);

return F;
}

void mode_solver::multiply_bloch_phase() {

std::vector<mpb_real> kvector = get_output_k();
Expand Down Expand Up @@ -1936,6 +2175,108 @@ std::vector<mpb_real> mode_solver::compute_group_velocity_component(vector3 d) {
return group_v;
}

/* as above, but only computes for given band */
mpb_real mode_solver::compute_1_group_velocity_component(vector3 d, int b) {
mpb_real u[3];
int ib = b - 1;
mpb_real group_v;
mpb_real scratch;

curfield_reset();

if (!mdata) {
meep::master_fprintf(stderr, "mode_solver.init must be called first!\n");
return group_v;
}

if (!kpoint_index) {
meep::master_fprintf(stderr, "mode_solver.solve_kpoint must be called first!\n");
return group_v;
}

/* convert d to unit vector in Cartesian coords: */
d = unit_vector3(matrix3x3_vector3_mult(Gm, d));
u[0] = d.x;
u[1] = d.y;
u[2] = d.z;

evectmatrix_resize(&W[0], 1, 0);
CHECK(nwork_alloc > 1, "eigensolver-nwork is too small");
evectmatrix_resize(&W[1], 1, 0);

maxwell_compute_H_from_B(mdata, H, W[1], (scalar_complex *) mdata->fft_data, ib, 0, 1);
maxwell_ucross_op(W[1], W[0], mdata, u);
evectmatrix_XtY_diag_real(W[1], W[0], &group_v, &scratch);

/* Reset scratch matrix sizes: */
evectmatrix_resize(&W[1], W[1].alloc_p, 0);
evectmatrix_resize(&W[0], W[0].alloc_p, 0);

if (freqs[ib] == 0) { /* v is undefined in this case */
group_v = 0.0; /* just set to zero */
}
else {
group_v /= negative_epsilon_ok ? sqrt(fabs(freqs[ib])) : freqs[ib];
}
return group_v;
}

/* returns group velocity for band b, in Cartesian coordinates */
vector3 mode_solver::compute_1_group_velocity(int b) {
vector3 v;
vector3 d;
matrix3x3 RmT = matrix3x3_transpose(Rm);
d.x = 1;
d.y = 0;
d.z = 0;
v.x = compute_1_group_velocity_component(matrix3x3_vector3_mult(RmT,d),b);
d.y = 1;
d.x = 0;
d.z = 0;
v.y = compute_1_group_velocity_component(matrix3x3_vector3_mult(RmT,d),b);
d.z = 1;
d.y = 0;
d.x = 0;
v.z = compute_1_group_velocity_component(matrix3x3_vector3_mult(RmT,d),b);

return v;
}

/* as above, but returns "group velocity" given by gradient of
frequency with respect to k in reciprocal coords ... this is useful
for band optimization. */
vector3 mode_solver::compute_1_group_velocity_reciprocal(int b) {
return matrix3x3_vector3_mult(matrix3x3_transpose(Gm),
compute_1_group_velocity(b));
}

/* compute the fraction of the field energy that is located in the
given range of dielectric constants: */
mpb_real mode_solver::compute_energy_in_dielectric(mpb_real eps_low, mpb_real eps_high) {
mpb_real *energy = (mpb_real *) curfield;
mpb_real epsilon = 0.0;
mpb_real energy_sum = 0.0;

if (!curfield || !strchr("DHBR", curfield_type)) {
meep::master_fprintf(stderr, "The D or H energy density must be loaded first.\n");
return 0.0;
}

int N = mdata->fft_output_size;

for (int i = 0; i < N; ++i) {
epsilon = mean_medium_from_matrix(mdata->eps_inv +i);
if (epsilon >= eps_low && epsilon <= eps_high) {
energy_sum += energy[i];
}
}
mpi_allreduce_1(&energy_sum, mpb_real, SCALAR_MPI_TYPE, MPI_SUM, mpb_comm);
energy_sum *= vol / H.N;
return energy_sum;
}



bool mode_solver::with_hermitian_epsilon() {
#ifdef WITH_HERMITIAN_EPSILON
return true;
Expand Down
Loading

0 comments on commit 4fdadb4

Please sign in to comment.