From 4fdadb436fc5fbdf3487744ee428bde349a676ec Mon Sep 17 00:00:00 2001 From: Christopher Hogan Date: Thu, 8 Mar 2018 20:33:46 -0600 Subject: [PATCH] Add missing MPB functionality (#223) * Add first_brillouin_zone and test * Add missing functionality * Add missing functionality --- libpympb/pympb.cpp | 361 +++++++++++++++++- libpympb/pympb.hpp | 20 + python/mpb.i | 87 +++++ python/solver.py | 126 +++++- python/tests/data/tutorial-d.k01.b02.te.h5 | Bin 0 -> 800 bytes python/tests/data/tutorial-d.k16.b02.te.h5 | Bin 0 -> 800 bytes python/tests/data/tutorial-mu.h5 | Bin 0 -> 10424 bytes .../data/tutorial-tot.rpwr.k16.b08.te.h5 | Bin 0 -> 10464 bytes python/tests/mpb.py | 134 ++++++- 9 files changed, 708 insertions(+), 20 deletions(-) create mode 100644 python/tests/data/tutorial-d.k01.b02.te.h5 create mode 100644 python/tests/data/tutorial-d.k16.b02.te.h5 create mode 100644 python/tests/data/tutorial-mu.h5 create mode 100644 python/tests/data/tutorial-tot.rpwr.k16.b08.te.h5 diff --git a/libpympb/pympb.cpp b/libpympb/pympb.cpp index b758b4102..e3d03b3de 100644 --- a/libpympb/pympb.cpp +++ b/libpympb/pympb.cpp @@ -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");\ @@ -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 @@ -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 = '-'; @@ -1595,6 +1672,168 @@ std::vector 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 kvector = get_output_k(); @@ -1936,6 +2175,108 @@ std::vector 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; diff --git a/libpympb/pympb.hpp b/libpympb/pympb.hpp index 97d966a39..bc93c3fe1 100644 --- a/libpympb/pympb.hpp +++ b/libpympb/pympb.hpp @@ -15,6 +15,7 @@ namespace py_mpb { struct mode_solver { static const int MAX_NWORK = 10; static const char epsilon_CURFIELD_TYPE = 'n'; + static const char mu_CURFIELD_TYPE = 'm'; static const int NUM_FFT_BANDS = 20; int num_bands; @@ -90,6 +91,7 @@ struct mode_solver { int get_kpoint_index(); void set_kpoint_index(int i); void get_epsilon(); + void get_mu(); void get_epsilon_tensor(int c1, int c2, int imag, int inv); void get_material_pt(meep_geom::material_type &material, vector3 p); void material_epsmu(meep_geom::material_type material, symmetric_matrix *epsmu, @@ -132,12 +134,30 @@ struct mode_solver { std::vector get_dims(); std::vector get_output_k(); + mpb_real 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); + mpb_real interp_val(vector3 p, int nx, int ny, int nz, int last_dim_size, + mpb_real *data, int stride, int conjugate); + scalar_complex interp_cval(vector3 p, int nx, int ny, int nz, int last_dim_size, + mpb_real *data, int stride); + symmetric_matrix interp_eps_inv(vector3 p); + + mpb_real get_epsilon_point(vector3 p); + cmatrix3x3 get_epsilon_inverse_tensor_point(vector3 p); + mpb_real get_energy_point(vector3 p); + cvector3 get_field_point(vector3 p); + cvector3 get_bloch_field_point(vector3 p); + void multiply_bloch_phase(); void fix_field_phase(); void compute_field_divergence(); std::vector compute_zparities(); std::vector compute_yparities(); std::vector compute_group_velocity_component(vector3 d); + mpb_real compute_1_group_velocity_component(vector3 d, int b); + vector3 compute_1_group_velocity(int b); + vector3 compute_1_group_velocity_reciprocal(int b); + mpb_real compute_energy_in_dielectric(mpb_real eps_low, mpb_real eps_high); bool with_hermitian_epsilon(); private: diff --git a/python/mpb.i b/python/mpb.i index 598626d35..2dd0a6a98 100644 --- a/python/mpb.i +++ b/python/mpb.i @@ -104,6 +104,65 @@ static int pylattice_to_lattice(PyObject *py_lat, lattice *l) { return 1; } + +static PyObject* v3_to_pyv3(vector3 *v) { + PyObject *geom_mod = PyImport_ImportModule("meep.geom"); + PyObject *v3_class = PyObject_GetAttrString(geom_mod, "Vector3"); + PyObject *args = Py_BuildValue("(ddd)", v->x, v->y, v->z); + PyObject *py_v = PyObject_Call(v3_class, args, NULL); + + Py_DECREF(geom_mod); + Py_DECREF(args); + Py_DECREF(v3_class); + + return py_v; +} + +static PyObject* cv3_to_pyv3(cvector3 *cv) { + PyObject *geom_mod = PyImport_ImportModule("meep.geom"); + PyObject *v3_class = PyObject_GetAttrString(geom_mod, "Vector3"); + + vector3 r = cvector3_re(*cv); + vector3 i = cvector3_im(*cv); + + Py_complex x, y, z; + x.real = r.x; + x.imag = i.x; + y.real = r.y; + y.imag = i.y; + z.real = r.z; + z.imag = i.z; + + PyObject *args = Py_BuildValue("(DDD)", &x, &y, &z); + PyObject *py_v = PyObject_Call(v3_class, args, NULL); + + Py_DECREF(geom_mod); + Py_DECREF(args); + Py_DECREF(v3_class); + + return py_v; +} + +static PyObject* cmatrix3x3_to_pymatrix(cmatrix3x3 *m) { + PyObject *c1 = cv3_to_pyv3(&m->c0); + PyObject *c2 = cv3_to_pyv3(&m->c1); + PyObject *c3 = cv3_to_pyv3(&m->c2); + + PyObject *geom_mod = PyImport_ImportModule("meep.geom"); + PyObject *matrix_class = PyObject_GetAttrString(geom_mod, "Matrix"); + + PyObject *args = Py_BuildValue("(OOO)", c1, c2, c3); + PyObject *res = PyObject_Call(matrix_class, args, NULL); + + Py_DECREF(c1); + Py_DECREF(c2); + Py_DECREF(c3); + Py_DECREF(geom_mod); + Py_DECREF(matrix_class); + Py_DECREF(args); + + return res; +} %} %include "std_string.i" @@ -132,6 +191,8 @@ static int pylattice_to_lattice(PyObject *py_lat, lattice *l) { meep_geom::material_data* }; +%apply double { mpb_real }; + %typemap(in) lattice { if (!pylattice_to_lattice($input, &$1)) { PyErr_PrintEx(0); @@ -169,6 +230,30 @@ static int pylattice_to_lattice(PyObject *py_lat, lattice *l) { } } +%typemap(out) vector3 { + $result = v3_to_pyv3(&$1); + + if (!$result) { + SWIG_fail; + } +} + +%typemap(out) cvector3 { + $result = cv3_to_pyv3(&$1); + + if (!$result) { + SWIG_fail; + } +} + +%typemap(out) cmatrix3x3 { + $result = cmatrix3x3_to_pymatrix(&$1); + + if (!$result) { + SWIG_fail; + } +} + %include "pympb.hpp" %pythoncode %{ @@ -193,6 +278,8 @@ static int pylattice_to_lattice(PyObject *py_lat, lattice *l) { output_charge_density, output_bpwr, output_dpwr, + output_tot_pwr, + output_dpwr_in_objects, output_poynting, output_poynting_x, output_poynting_y, diff --git a/python/solver.py b/python/solver.py index 588796ae4..7575db258 100644 --- a/python/solver.py +++ b/python/solver.py @@ -1,5 +1,6 @@ from __future__ import division, print_function +import functools import os import numbers import re @@ -146,6 +147,9 @@ def ExH(e, h): return flat_res + def get_epsilon(self): + self.mode_solver.get_epsilon() + def get_bfield(self, which_band): return self._get_field('b', which_band) @@ -180,6 +184,44 @@ def _get_field(self, f, band): return field + def get_epsilon_point(self, p): + return self.mode_solver.get_epsilon_point(p) + + def get_epsilon_inverse_tensor_point(self, p): + return self.mode_solver.get_epsilon_inverse_tensor_point(p) + + def get_energy_point(self, p): + return self.mode_solver.get_energy_point(p) + + def get_field_point(self, p): + return self.mode_solver.get_field_point(p) + + def get_bloch_field_point(self, p): + return self.mode_solver.get_bloch_field_point(p) + + def get_tot_pwr(self, which_band): + self.get_dfield(which_band) + self.compute_field_energy() + + dims = self.mode_solver.get_dims() + epwr = np.zeros(np.prod(dims)) + self.mode_solver.get_curfield(epwr) + epwr = np.reshape(epwr, dims) + + self.get_bfield(which_band) + self.compute_field_energy() + + hpwr = np.zeros(np.prod(dims)) + self.mode_solver.get_curfield(hpwr) + hpwr = np.reshape(hpwr, dims) + + tot_pwr = epwr + hpwr + + self.mode_solver.set_curfield(tot_pwr.ravel()) + self.mode_solver.set_curfield_type('R') + + return tot_pwr + # The band-range-data is a list of tuples, each consisting of a (min, k-point) # tuple and a (max, k-point) tuple, with each min/max pair describing the # frequency range of a band and the k-points where it achieves its minimum/maximum. @@ -288,10 +330,8 @@ def output_epsilon(self): self.output_field_to_file(mp.ALL, self.get_filename_prefix()) def output_mu(self): - print("output_mu: Not yet supported") - # TODO - # self.mode_solver.get_mu() - # self.mode_solver.output_field_to_file(-1, self.get_filename_prefix) + self.mode_solver.get_mu() + self.output_field_to_file(mp.ALL, self.get_filename_prefix()) def output_field_to_file(self, component, fname_prefix): curfield_type = self.mode_solver.get_curfield_type() @@ -445,9 +485,15 @@ def _create_fname(self, fname, prefix, parity_suffix): def compute_field_energy(self): return self.mode_solver.compute_field_energy() + def compute_field_divergence(self): + mode_solver.compute_field_divergence() + def compute_energy_in_objects(self, objs): return self.mode_solver.compute_energy_in_objects(objs) + def compute_energy_in_dielectric(self, eps_low, eps_high): + return self.mode_solver.compute_energy_in_dielectric(eps_low, eps_high) + def compute_group_velocities(self): xarg = mp.cartesian_to_reciprocal(mp.Vector3(1), self.geometry_lattice) vx = self.mode_solver.compute_group_velocity_component(xarg) @@ -458,6 +504,16 @@ def compute_group_velocities(self): return [mp.Vector3(x, y, z) for x, y, z in zip(vx, vy, vz)] + def compute_group_velocity_component(self, direction): + return self.mode_solver.compute_group_velocity_component(direction) + + def compute_one_group_velocity(self, which_band): + return self.mode_solver.compute_1_group_velocity(which_band) + + def compute_one_group_velocity_component(self, direction, which_band): + return self.mode_solver.compute_1_group_velocity_component(direction, + which_band) + def randomize_fields(self): self.mode_solver.randomize_fields() @@ -715,6 +771,42 @@ def bfunc(ms, b_prime): return ks + def first_brillouin_zone(self, k): + """ + Function to convert a k-point k into an equivalent point in the + first Brillouin zone (not necessarily the irreducible Brillouin zone) + """ + def n(k): + return mp.reciprocal_to_cartesian(k, self.geometry_lattice).norm() + + def try_plus(k, v): + if n(k + v) < n(k): + return try_plus(k + v, v) + else: + return k + + def _try(k, v): + return try_plus(try_plus(k, v), mp.Vector3() - v) + + try_list = [ + mp.Vector3(1, 0, 0), mp.Vector3(0, 1, 0), mp.Vector3(0, 0, 1), + mp.Vector3(0, 1, 1), mp.Vector3(1, 0, 1), mp.Vector3(1, 1, 0), + mp.Vector3(0, 1, -1), mp.Vector3(1, 0, -1), mp.Vector3(1, -1, 0), + mp.Vector3(1, 1, 1), mp.Vector3(-1, 1, 1), mp.Vector3(1, -1, 1), + mp.Vector3(1, 1, -1), + ] + + def try_all(k): + return functools.reduce(_try, try_list, k) + + def try_all_and_repeat(k): + knew = try_all(k) + return try_all_and_repeat(knew) if n(knew) < n(k) else k + + k0 = k - mp.Vector3(*[round(x) for x in k]) + + return try_all_and_repeat(k0) if n(k0) < n(k) else try_all_and_repeat(k) + # Predefined output functions (functions of the band index), for passing to `run` @@ -810,6 +902,32 @@ def output_dpwr(ms, which_band): ms.output_field() +def output_tot_pwr(ms, which_band): + ms.get_tot_pwr(which_band) + ms.output_field_to_file(-1, ms.get_filename_prefix() + 'tot.') + + +def output_dpwr_in_objects(output_func, min_energy, objects=[]): + """ + The following function returns an output function that calls output_func for + bands with D energy in objects > min-energy. For example, + output_dpwr_in_objects(output_dfield, 0.20, some_object) would return an + output function that would spit out the D field for bands with at least %20 + of their D energy in some-object. + """ + + def _output(ms, which_band): + ms.get_dfield(which_band) + ms.compute_field_energy() + energy = ms.compute_energy_in_objects(objects) + fmt = "dpwr:, {}, {}, {} " + print(fmt.format(which_band, ms.freqs[which_band - 1], energy)) + if energy >= min_energy: + apply_band_func(ms, output_func, which_band) + + return _output + + def output_charge_density(ms, which_band): ms.get_charge_density(which_band) ms.output_field_to_file(-1, ms.get_filename_prefix()) diff --git a/python/tests/data/tutorial-d.k01.b02.te.h5 b/python/tests/data/tutorial-d.k01.b02.te.h5 new file mode 100644 index 0000000000000000000000000000000000000000..82daee09f5d9bab4d4ef6cf1c24d6067e7457f8a GIT binary patch literal 800 zcmeD5aB<`1lHy_j0S*oZ76t(@6Gr@p0tIG>2#gPtPk=HQp>zk7Ucm%mFfxE31A_!q zTo7tLy1I}cS67e{nE5aos%?}S;UVDR>KFhDf(U3ha6su3&~ygng3}s^4OR>jq<{th Dq@*Z_ literal 0 HcmV?d00001 diff --git a/python/tests/data/tutorial-d.k16.b02.te.h5 b/python/tests/data/tutorial-d.k16.b02.te.h5 new file mode 100644 index 0000000000000000000000000000000000000000..82daee09f5d9bab4d4ef6cf1c24d6067e7457f8a GIT binary patch literal 800 zcmeD5aB<`1lHy_j0S*oZ76t(@6Gr@p0tIG>2#gPtPk=HQp>zk7Ucm%mFfxE31A_!q zTo7tLy1I}cS67e{nE5aos%?}S;UVDR>KFhDf(U3ha6su3&~ygng3}s^4OR>jq<{th Dq@*Z_ literal 0 HcmV?d00001 diff --git a/python/tests/data/tutorial-mu.h5 b/python/tests/data/tutorial-mu.h5 new file mode 100644 index 0000000000000000000000000000000000000000..bbba29cdf88f0b83645e8c739b617317888be23a GIT binary patch literal 10424 zcmeHJL2DC16n?u&>k<)6qM)8udr%OhhZa3Vc7;Yl6&3X$Vxeog;6mGkbd~DCUgJOT zqzHmNdJ&5U6``e2yb9j@18TsNUg}KVJaMwb4M|fs;cdvg*?Dii%=hM-mw8m0I=?S- zBm?rXEJ(xtW|x|6^{k!15qc;8tm+ub)|Fi~v15Sz031@|S%u#kcl~Oy2qKd3$*O5E zp2nWOI`Xm|_6(GY=dMVOHx$uvWjB@1)GN$+p(p$EUKj>t-@fOU!%A%l=KQ5{EvSY; zWw9AbdR4Nf=@DzdXj8w~NQL4I!$bsI&1h#+LisW>Ydc^yV_lfo3lrNZoH!yzrQb9W z&{FX!V-~{!Y>dPdaYWnF-fIEM?``*8H>GG5TDYXi*IK`vg58KL7A8i+b<8-c< zNB$i*EZ#U>QCi#Eu62b1hLV=b|E;TNK2ZKr=kk$?ll{E-O0Bd%H>-*tDM|#3C(7?J zI=J113^W}45w@HMnE}=T>wtB@=YZ<~*8#2rTnBou17n|V{rGz3mt)jEK6!I7`PGHopLg<8{Ggt{%JG#CXAd@Xn|&z_a2#6SlBqBx}$(RrovxtHi zNFqToI|!&`Mi7BZ21WdRhTi}E_owE)s$cc$SG{^zdhg+!vv>F2y}MVh<$m)8X3`QW z5*+qNOpGJMk^AE-{pUOW)M||UhF;VA_V|U1pON^PCNT1a%VGa;ILi2Yc^vWWPux({mjt{Ewgece(q| z;)*WoUw{6e<7y;;Ip%-2l>gsx_2=~%34b=u*&jCI*zXGPcQ*X}>vNNa{;^N}?dQ+_ z#sBT*f4=?8RzRSr?6l`Ea?HWERiV<$<(M%yB^F19$TGQI++C-#r5XL1a|ZX>OEPVo zGnW>|i7^GYgNFBn3p06HhgUdR2r$8Yvkq1n4HK2w4hMHl7$AGpvf2txc99~b@q0AG zTS=5`eU@W=J^4!JZRosRNnDjKye^;jlAQ5(SS5M0h@{3W7kLw1Kr-;!8}F+yBhMee z_X}{m5FCFr&bJHa&p~@0qJ8aXuL0Um{TYq^QGXRq-}u9Cm*?w`?wBaYw25De6H%9C zPRer}?m9~|uJ*Fb&M}fqe|=SYyqy?h7E@w$-BOrYEdI*PQcHlbPSP9nRvacxv&cA? z$pfU#Y@u6Hd>0`*)t8FKw~}hVCq60@8j0cAXXhvWUPUy!L;P(tUy{XNwT(=zN=VeQ zvZx@9M?~2BtH>p|jxny=I128IM7te$*;GZySuX!Kcwd0`#yTmkl|1YlGD(2(6`$I_ zwSoD9%|k;%VgC4Nsh*`^&-=`QhD@+;c=h+A$H3lI_m<@cgZ(WNMOSBnKQ+0HKH}hC z#^;90K=Ai_@^lM!zt7)X6wmq{xOKWs-*_2DmLYIrViBSA%t->jm0~*<6y% z5pE*Oy>}|6nbr`+L33lH=Wj`M>=dUHj<1MpoI-kH_zRLQw57-du7#IW{h9~&JG8#I z&o!$drEWWpCkQr?eCgNA37BNPI`BWmZUxeA6F4o%HBFG3IIvydY3;tUd*ZU6*5iQG7 z4TH9R;=T9>{%!qkV!5k2V0c4238~TBtm@uOG%984cIT`exp0aDVsJVR0XLe(Lg*+WH>NWS~^g=H2>svX}R%_8N?T zYP35~8s_^j)j;<_W z*l#|{uEP-rn@wyy411eTiG=graz)>cjdSrv?OJxf3!Z8#uznwJ36nLQCd%ADe9~;8 ztS}?MtKIRTP=N7xU>0)6fx~zSJGd0943Yj`>AIEY`^mc)X~SDDy9wiLXkfdIPtH3< z-x50giSU~C8b?TdB#H{lc5tsZl7OZ82X3FQC&%79q`n`hBh@0-woY(uezkgqKit1B zxqkCx>5s(odVuZj6Q4-Uk~7K5ethDc`lbIFj9=)r`Oq<#@1en(%ovz|_Cx2w3b5zJ z1*OhMVBh5&=XQIrcd@n69bK@$Voub2b{zx?Iv%m{FgIc8JN8_R&vxl#^IGF|X~_gQ zFMOGxAvVtCLe;sVu-|HB@7J?_%S!sC&RHeMC~&U+^vQ+GJe{t8nO{3h+8<5HkNWzB z1kZcVbbI%c+x*YXqm_HeqH8ZThFv`U0 zD*Z_ATBT*j#AmCRJIcOXbdcoWzk^hUCc zjfaDJwLi1x;;!xPPi%g$E-}gY1@gJB`uI{duf3XUaE?9SLD_lceGum_Dx?b8{qC?I zwV3rgqU>2(XT^8o7^ELumj`jLaLuRw$slp$MKqo1?jzTgo~(7v?;)>y?tQ*l)`-i@J&@VSJm+ppC@athjVwd&f5R%KDqDY*V6TK=dBk52@i zC=JBH^VZg18mjX7M0exfZShaLh_#m3sS7YZ|G`W5L6|R8(Pr5-nBVoB+*+gjA)${A5`2_2?qRer#11G)^i6<>H*ZK|+ z)w@g7C!FmgCa0g&Ro?F**B14r%kJ+c^L?7G{OHk5>L({<2$*z{hzElOY3Cp=q7okH zHFS~*ueZhO^E=6v`YIoHxZbp0MSUmS&$=dIcXUn{k^A|qo2XYe(J8O-@rvvwy6!nA zb71`bSw^YDDSc#3d_-6-%)i2Xe3KQ}quo&K&#r^8ORvsdus3|sLW4@M-+JJ&ARqFB z(!ERdY@U=@70t5%e~(!<_*6sv{(R`Pkss6%Id|KSS&&ctemtQb2zf=h@PS4toUbE+ z$$Jk#oRw7B$Fln^tWqz?`Yjlfo_}gyKk*tX-LBc%OV%cC8}Dz^OLjgp-4YenLwK{u zr}+!JiO0)z!=C$k$nhvZTVmV`Lb#YdbiA$mb= z(Mfy2ZkOJ^yfNo{i1g%($ub>1#BA?xzv|uZB@bTbKCrCmCaTS|rZZ_hWbE(vEDZhn zNPyELF|o!zqU1L`#vQKnWUL;{fcv!vH_WZA?k1621%y;{ddY$X_Q@^nJ>-kKi1piZ zJtX&jm~|Y?CscJVcR#xha>7f#K|IL2%9gX|qCfZe#6qxlfzTT97_firg`;ES!5@CS z=(Gs%&)X?u-+J&@=yBTB$AH(oY(4j<0Is-lKjVrA)YtIyG8*@w&dhnSSpE*=TlMLA zWn&=k#0#5vv*)}0w6akJ#Gz}&D>Gl%kF23;Cq?l4N^MA4aepu2UFu6NU)4jZe(g-C zdEGt}NwY%YfT6{qg) z+m!O196V$__Yk}gq+DRte(BHx4wO+?9W-@TT-3*2F!y(Q? zTJM>N!+v|$Ule!-c2=wv&atcSB`Tdkk1Nk}8O^?RP71}MjQWu~OEa^@nS^)wZK~#y z%=RC|%GM7^G6o(UO?GhYI<@iBLO8Q`;(+}&5K;7T{B^YXhSavj!Zjp z3dUa%`P{z|;$gqtrTUd{E=r$2W;`K3IC_gdii13<9bUd~KGX~2V}lVBpl;ec%YWYj z_~EwN<6$$vQQ-xO9ajNQzTKPJwiEDqowM9CJ;0f}x*{A$1753hQO(i;+~RaYO?4F1 z@A5#y%w15&@^ee0Jt3bPBqq(90eP)3p=pB&obN1Q85IqPbKkV6iYKsN9Zp>HV(@!^ zMSDZWLUD#uGk2TZc5%k*#80{Mt`f}4ysk5y^3u#E|2M%yPitnLbEXaJQoD6ldCdr{P3h={Tu@K*dfTN;JYL;FDJSlYWjA;k>>vXiwd}?MT zxg#&Gv=Q*R>GHJIzX49Iw3GZc1o&%b(_PtL0k@RLthu=p>bIV2;@b&Omn?7E3kX3z zFV5R>!5;Ej-oVO}tKfVm&oloT3vsTS@xt>7>|cjQ(@z}%;*9ybD`(vjB$(yRqK=!+ z%Q4=YT~aL$C^3aKL;S^!%8bn}EgLDg{=6(>at+)M^IC8uHc^6^-{$c~4A$Y@>8Iwg zFy8H{uA#zkF0y~r)2@d6@cXpp!{U%9%Y2J7qM=?etrwhQ!9I`lhgPM4AA&N2Reu5; zW$nAw;Q;ufwZ`hO7vO4L0J8Gn@2pSGM@G+oNmNVDu1ti!qbp?ve;eR$l_`V%a!Wsu zqPELDih#!+6^Bm^1l;o0Dpl$`)bA_buSi^gI%RfZ?OF%OLn6weQ}iKkd2Ah>IR(z| zmiO@*b0E%dSA`80x=1klTD1HYgI_t47q!nijc2aL96S16U5%OMBR*&18g<4&K&(|= zLyZw?@GTv19MAmT=&4--&nv3sguH=x>0e=1coD|aou=pJ3FD`1>`qjL`AYZ=f=6Kf z@NV;sJHVdhF5E!|@Pn&D+ACqOck$#W4Uu4fc+t9jwt%Z=kBi?t2>$u`sQqXK{`QL= zwKxQQXo-pKh?C+M8 zZ`S4Dx9{!a2cqBo$OwuoS-!|xhuL$k^{G8v?-onWoCNFebFJx5O;Jirig@XTdGLO4 zi%D4sjMu$#>fw7Z{?aauA`h6aFl+VwqcHymDW3vQuqWncS9KOg-E5kd;tuxePE!xv z4)#Bjnaz6+{;1@sB?W?i$(FB0O{`1EIKIugjmHW|vBCLc($gQ2Q2~7Wa_Gl4s68Fm z4t+_&#h~<7z~{TH^omUYr>bd<<4y&<_ExCv5CPnh`26~frBJ`6nlBovLmeAi`mMqq z^7))-{!_hvN#@3hf)BryN;5*`GtNwbIM?(~$Wws*F1oziVk!9DQrg|pt)l76C6(k$Ct(V=={{4>1f{`v^NIrr~XKwe?{o;SKR;exE>R6eO}>u>EQa^!1bi- zTZQXQ*WUv7qZjul2>0s@?%!XAi=Z&5}dLHTd?8Wm+ z&us>k4e_9U@ zVSS|avJUI#e5|L>u)flI>x=dG6xQP;tk30Gul=!p3t>G^!1~^R^`6%MSBM9E5g!O5 zUWh>aK=Fhh;)`2|H;NE{=p!B}Kz#BW;uQzPFExl~f)L+SBi?C5{1<_EC>!yS9pa@q zh@U8)qWG!^@zz1aUlfl~e8xw-7Kr#w7xA1s;yZoBdx41mk`NDSAU>pcaU9~uD#Vk1 zh%YJLbVd9*fOvEY;!}!O;}E}oMLf$OzEwfIOYwgl;^8NVk82PwYao86csdO6bpzsU zTg2Zfh{q{Dr+7UQ@p~2G`4@=qe?`1M1@ZqJ><1jMKZwVEf%XrypSXnmMJDze&e(sn zVLuXw{R!%|_Ro{CpVqc*!HOQCjMgBw*`4nB`S1gcku|xi42>BSw&rrT*6!JIO z$mcj9zjF=wo{7l+3?d&CiTn`di{>GJMEN8J`6bFXNg@9<7Wt?Z$WKweYBBOx`;pIT zM1G6%U6lX3g?t$0$3l@Wi$ng*9r-jpArs zzK`;MpOFvzi2Ps}@`b9%AHG38@domX4#+oB{?Qrv$Z+H*`N&sN{*v;Ul;5O$=Q!kl zbC3@`fc)rAIt=|FZ7|_K=lVN)FXsYpLl?J1=TMKQO^)VeS_*9RR1}OddMu) zM|4pyS%>-w)l>3NU!i(S9qKP9QI82jeWnlf8V%HM{7}z1j{43W)O!k1{}KIRL=U=) z`j9^AMSD;`T7`O&FX~I0s5jZ5{zUbtJ*ZC!p(A?eXaoYx+2u?zM-D?`{EIO?;PrV(y0G^L_JUk z^}(M}FEmH}a0}{*!l*Agqu#g<^~ckwM^b$<67@=D)Gq~5&-6!qlj@zLQU9fSXdUXK z&8U}B{q!s9snw{j#-iRDf%>Zw>ahi=&kCYmYmEBs*P|nP?qSq-2e*#sy^*N@TB07T zkNPk{y|^CrV`J2l+fiSxLA}`u_2(t1M^k%d2Ms^+o8$Z;=R?)%a5lsQ*O;> NM@!pKceXC)zW_)l