Skip to content

Commit

Permalink
Fix performance issues in ModeSolver.find_k (NanoComp#644)
Browse files Browse the repository at this point in the history 
* Reuse fields if able

* Allow reusing ModeSolver objects

* Don't prematurely evaluate rootfun

* Add property accessor for ModeSolver's C++ mode_solver

* Fix resolution and lattice accessors

* Adjust grid_size when resolution or geometry_lattice are changed
  • Loading branch information
@ChristopherHogan @stevengj
ChristopherHogan authored and stevengj committed Dec 31, 2018
1 parent cf17fae commit caef3c5
Show file tree
Hide file tree
Showing 5 changed files with 366 additions and 108 deletions.
122 changes: 91 additions & 31 deletions libpympb/pympb.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -185,14 +185,11 @@ static void deflation_constraint(evectmatrix X, void *data) {
/******* mode_solver *******/

mode_solver::mode_solver(int num_bands,
int parity,
double resolution[3],
lattice lat,
double tolerance,
int mesh_size,
meep_geom::material_data *_default_material,
geometric_object_list geom,
bool reset_fields,
bool deterministic,
double target_freq,
int dims,
Expand All @@ -208,7 +205,6 @@ mode_solver::mode_solver(int num_bands,
int eigensolver_nwork,
int eigensolver_block_size):
num_bands(num_bands),
parity(parity),
target_freq(target_freq),
tolerance(tolerance),
mesh_size(mesh_size),
Expand All @@ -218,6 +214,7 @@ mode_solver::mode_solver(int num_bands,
force_mu(force_mu),
use_simple_preconditioner(use_simple_preconditioner),
grid_size(grid_size),
nwork_alloc(0),
eigensolver_nwork(eigensolver_nwork),
eigensolver_block_size(eigensolver_block_size),
last_parity(-2),
Expand All @@ -235,11 +232,17 @@ mode_solver::mode_solver(int num_bands,
curfield_type('-'),
eps(true) {

this->lat = lat;

geometry_lattice = lat;
dimensions = dims;
ensure_periodicity = periodicity;
geometry_tree = NULL;
H.data = NULL;
Hblock.data = NULL;
muinvH.data = NULL;

for (int i = 0; i < MAX_NWORK; i++) {
W[i].data = NULL;
}

for (int i = 0; i < 3; ++i) {
this->resolution[i] = resolution[i];
Expand All @@ -249,21 +252,14 @@ mode_solver::mode_solver(int num_bands,
}
}

default_material = _default_material;

#ifndef WITH_HERMITIAN_EPSILON
meep_geom::medium_struct *m;
if (meep_geom::is_medium(default_material, &m)) {
if (meep_geom::is_medium(_default_material, &m)) {
meep_geom::check_offdiag(m);
}
#else
(void)_default_material;
#endif

geometry = geom;

// `init` is called in the constructor to avoid the need to copy the
// geometric objects and default material. They can then be safely freed by
// typemaps once the mode_solver constructor call returns to python.
init(parity, reset_fields);
}

mode_solver::~mode_solver() {
Expand Down Expand Up @@ -716,9 +712,12 @@ bool mode_solver::using_mu() {
return mdata && mdata->mu_inv != NULL;
}

void mode_solver::init(int p, bool reset_fields) {
void mode_solver::init(int p, bool reset_fields, geometric_object_list geometry,
meep_geom::material_data *_default_material) {
int have_old_fields = 0;

default_material = _default_material;

n[0] = grid_size.x;
n[1] = grid_size.y;
n[2] = grid_size.z;
Expand Down Expand Up @@ -759,6 +758,36 @@ void mode_solver::init(int p, bool reset_fields) {
block_size = num_bands;
}

if (mdata) {
if (n[0] == mdata->nx && n[1] == mdata->ny && n[2] == mdata->nz &&
block_size == Hblock.alloc_p && num_bands == H.p &&
eigensolver_nwork + (mdata->mu_inv != NULL) == nwork_alloc) {

have_old_fields = 1;
}
else {
destroy_evectmatrix(H);
for (int i = 0; i < nwork_alloc; ++i) {
destroy_evectmatrix(W[i]);
}

if (Hblock.data != H.data) {
destroy_evectmatrix(Hblock);
}
if (muinvH.data != H.data) {
destroy_evectmatrix(muinvH);
}
}
destroy_maxwell_target_data(mtdata);
mtdata = NULL;
destroy_maxwell_data(mdata);
mdata = NULL;
curfield_reset();
}
else {
srand(time(NULL));
}

if (deterministic) {
// seed should be the same for each run, although
// it should be different for each process.
Expand All @@ -774,7 +803,7 @@ void mode_solver::init(int p, bool reset_fields) {
mtdata = create_maxwell_target_data(mdata, target_freq);
}

init_epsilon();
init_epsilon(&geometry);

if (check_maxwell_dielectric(mdata, 0)) {
meep::abort("invalid dielectric function for MPB");
Expand Down Expand Up @@ -816,7 +845,7 @@ void mode_solver::init(int p, bool reset_fields) {
evectmatrix_flops = eigensolver_flops;
}

void mode_solver::init_epsilon() {
void mode_solver::init_epsilon(geometric_object_list *geometry) {
mpb_real no_size_x = geometry_lattice.size.x == 0 ? 1 : geometry_lattice.size.x;
mpb_real no_size_y = geometry_lattice.size.y == 0 ? 1 : geometry_lattice.size.y;
mpb_real no_size_z = geometry_lattice.size.z == 0 ? 1 : geometry_lattice.size.z;
Expand Down Expand Up @@ -844,20 +873,20 @@ void mode_solver::init_epsilon() {
matrix3x3_to_arr(R, Rm);
matrix3x3_to_arr(G, Gm);

geom_fix_objects0(geometry);
geom_fix_objects0(*geometry);

meep::master_printf("Geometric objects:\n");
if (meep::am_master()) {
for (int i = 0; i < geometry.num_items; ++i) {
for (int i = 0; i < geometry->num_items; ++i) {

#ifndef WITH_HERMITIAN_EPSILON
meep_geom::medium_struct *mm;
if (meep_geom::is_medium(geometry.items[i].material, &mm)) {
if (meep_geom::is_medium(geometry->items[i].material, &mm)) {
meep_geom::check_offdiag(mm);
}
#endif

display_geometric_object_info(5, geometry.items[i]);
display_geometric_object_info(5, geometry->items[i]);

// meep_geom::medium_struct *mm;
// if (meep_geom::is_medium(geometry.items[i].material, &mm)) {
Expand All @@ -884,7 +913,7 @@ void mode_solver::init_epsilon() {
b0.high.x += tmp_size.x / mdata->nx;
b0.high.y += tmp_size.y / mdata->ny;
b0.high.z += tmp_size.z / mdata->nz;
geometry_tree = create_geom_box_tree0(geometry, b0);
geometry_tree = create_geom_box_tree0(*geometry, b0);
}

if (verbose && meep::am_master()) {
Expand All @@ -896,14 +925,14 @@ void mode_solver::init_epsilon() {
int tree_nobjects;
geom_box_tree_stats(geometry_tree, &tree_depth, &tree_nobjects);
meep::master_printf("Geometric object tree has depth %d and %d object nodes"
" (vs. %d actual objects)\n", tree_depth, tree_nobjects, geometry.num_items);
" (vs. %d actual objects)\n", tree_depth, tree_nobjects, geometry->num_items);

// restricted_tree = geometry_tree;

reset_epsilon();
reset_epsilon(geometry);
}

void mode_solver::reset_epsilon() {
void mode_solver::reset_epsilon(geometric_object_list *geometry) {
int mesh[3] = {
mesh_size,
(dimensions > 1) ? mesh_size : 1,
Expand All @@ -922,7 +951,7 @@ void mode_solver::reset_epsilon() {
set_maxwell_dielectric(mdata, mesh, R, G, dielectric_function, mean_epsilon_func,
static_cast<void *>(this));

if (has_mu()) {
if (has_mu(geometry)) {
meep::master_printf("Initializing mu function...\n");
eps = false;
set_maxwell_mu(mdata, mesh, R, G, dielectric_function, mean_epsilon_func,
Expand All @@ -931,14 +960,14 @@ void mode_solver::reset_epsilon() {
}
}

bool mode_solver::has_mu() {
bool mode_solver::has_mu(geometric_object_list *geometry) {
// TODO: mu_file_func
if (material_has_mu(default_material) || force_mu) {
return true;
}

for (int i = 0; i < geometry.num_items; ++i) {
if (material_has_mu(geometry.items[i].material)) {
for (int i = 0; i < geometry->num_items; ++i) {
if (material_has_mu(geometry->items[i].material)) {
return true;
}
}
Expand Down Expand Up @@ -1001,6 +1030,11 @@ void mode_solver::set_parity(integer p) {
set_kpoint_index(0); /* reset index */
}

void mode_solver::set_num_bands(int nb) {
num_bands = nb;
freqs.resize(nb);
}

int mode_solver::get_kpoint_index() {
return kpoint_index;
}
Expand Down Expand Up @@ -1586,6 +1620,32 @@ std::vector<int> mode_solver::get_dims() {
return dims;
}

void mode_solver::set_grid_size(vector3 gs) {
grid_size.x = gs.x;
grid_size.y = gs.y;
grid_size.z = gs.z;
}

int mode_solver::get_libctl_dimensions() {
return dimensions;
}

void mode_solver::set_libctl_dimensions(int val) {
dimensions = val;
}

bool mode_solver::get_libctl_ensure_periodicity() {
return ensure_periodicity;
}

void mode_solver::set_libctl_ensure_periodicity(bool val) {
ensure_periodicity = val;
}

void mode_solver::set_libctl_geometry_lattice(lattice val) {
geometry_lattice = val;
}

void mode_solver::get_curfield(double *data, int size) {
mpb_real *p = (mpb_real *)curfield;

Expand Down
48 changes: 32 additions & 16 deletions libpympb/pympb.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,6 @@ struct mode_solver {
static const int NUM_FFT_BANDS = 20;

int num_bands;
int parity;
double resolution[3];
double target_freq;
lattice lat;
Expand Down Expand Up @@ -57,10 +56,7 @@ struct mode_solver {
int iterations;
double eigensolver_flops;

geometric_object_list geometry;

geom_box_tree geometry_tree;
// geom_box_tree restricted_tree;

mpb_real vol;

Expand All @@ -81,25 +77,39 @@ struct mode_solver {
evectmatrix muinvH;
evectmatrix W[MAX_NWORK];

meep_geom::material_data *default_md;

std::vector<mpb_real> freqs;

bool verbose;
bool deterministic;

mode_solver(int num_bands, int parity, double resolution[3], lattice lat, double tolerance,
int mesh_size, meep_geom::material_data *_default_material, geometric_object_list geom,
bool reset_fields, bool deterministic, double target_freq, int dims, bool verbose,
bool periodicity, double flops, bool negative_epsilon_ok, std::string epsilon_input_file,
std::string mu_input_file, bool force_mu, bool use_simple_preconditioner, vector3 grid_size,
int eigensolver_nwork, int eigensolver_block_size);
mode_solver(int num_bands,
double resolution[3],
lattice lat,
double tolerance,
int mesh_size,
meep_geom::material_data *_default_material,
bool deterministic,
double target_freq,
int dims,
bool verbose,
bool periodicity,
double flops,
bool negative_epsilon_ok,
std::string epsilon_input_file,
std::string mu_input_file,
bool force_mu,
bool use_simple_preconditioner,
vector3 grid_size,
int eigensolver_nwork,
int eigensolver_block_size);

~mode_solver();

void init(int p, bool reset_fields);
void init(int p, bool reset_fields, geometric_object_list geometry, meep_geom::material_data *_default_material);
void solve_kpoint(vector3 kpoint);
bool using_mu();
void set_parity(int p);
void set_num_bands(int nb);
int get_kpoint_index();
void set_kpoint_index(int i);
void get_epsilon();
Expand All @@ -112,9 +122,9 @@ struct mode_solver {
mpb_real d1, mpb_real d2, mpb_real d3, mpb_real tol, const mpb_real r[3]);

void randomize_fields();
void init_epsilon();
void reset_epsilon();
bool has_mu();
void init_epsilon(geometric_object_list *geometry);
void reset_epsilon(geometric_object_list *geometry);
bool has_mu(geometric_object_list *geometry);
bool material_has_mu(void *mt);
void curfield_reset();

Expand Down Expand Up @@ -146,6 +156,12 @@ struct mode_solver {
void set_curfield_type(char t);
std::string get_parity_string();
std::vector<int> get_dims();
void set_grid_size(vector3 gs);
int get_libctl_dimensions();
void set_libctl_dimensions(int val);
bool get_libctl_ensure_periodicity();
void set_libctl_ensure_periodicity(bool val);
void set_libctl_geometry_lattice(lattice val);
std::vector<mpb_real> get_output_k();

mpb_real get_val(int ix, int iy, int iz, int nx, int ny, int nz, int last_dim_size,
Expand Down
3 changes: 1 addition & 2 deletions python/geom.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -752,8 +752,7 @@ def newton(x, a, b, dx):
a_prime = x if f < 0 else a
b_prime = x if f > 0 else b

cond = f_memo(lazy(a_prime))[0] * f_memo(lazy(b_prime))[0] > 0
if dx != x_max - x_min and dx * (f / df) < 0 and cond:
if dx != x_max - x_min and dx * (f / df) < 0 and f_memo(lazy(a_prime))[0] * f_memo(lazy(b_prime))[0] > 0:
raise ValueError("failed to bracket the root in find_root_deriv")

if isinstance(a, numbers.Number) and isinstance(b, numbers.Number):
Expand Down
Loading

0 comments on commit caef3c5

Please sign in to comment.