Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix performance issues in ModeSolver.find_k #644

Merged
merged 6 commits into from
Dec 31, 2018
Merged

Fix performance issues in ModeSolver.find_k #644

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
3e2785b
Reuse fields if able
ChristopherHogan Dec 19, 2018
d59992f
Allow reusing ModeSolver objects
ChristopherHogan Dec 20, 2018
4363947
Don't prematurely evaluate rootfun
ChristopherHogan Dec 21, 2018
fd6f318
Add property accessor for ModeSolver's C++ mode_solver
ChristopherHogan Dec 29, 2018
7c303e7
Fix resolution and lattice accessors
ChristopherHogan Dec 31, 2018
06c88ab
Adjust grid_size when resolution or geometry_lattice are changed
ChristopherHogan Dec 31, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
122 changes: 91 additions & 31 deletions libpympb/pympb.cpp
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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