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

Jump to bottom

fix segmentation fault in material_grids_addgradient of adjoint solver #1574

Merged
merged 1 commit into from
May 19, 2021
Merged

fix segmentation fault in material_grids_addgradient of adjoint solver #1574

Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
34 changes: 17 additions & 17 deletions src/meepgeom.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -420,6 +420,14 @@ double material_grid_val(vector3 p, material_data *md) {

}

static double tanh_projection(double u, double beta, double eta) {
double u_proj = 0;
double tanh_beta_eta = tanh(beta*eta);
u_proj = (tanh_beta_eta + tanh(beta*(u-eta))) /
(tanh_beta_eta + tanh(beta*(1-eta)));
return u_proj;
}

double matgrid_val(vector3 p, geom_box_tree tp, int oi, material_data *md) {
double uprod = 1.0, umin = 1.0, usum = 0.0, udefault = 0.0, u;
int matgrid_val_count = 0;
Expand Down Expand Up @@ -458,15 +466,7 @@ double matgrid_val(vector3 p, geom_box_tree tp, int oi, material_data *md) {
: (md->material_grid_kinds == material_data::U_MEAN ? usum / matgrid_val_count
: udefault)));

// project interpolated grid point
double u_proj;
if (md->beta != 0) {
double tanh_beta_eta = tanh(md->beta*md->eta);
u_proj = (tanh_beta_eta + tanh(md->beta*(u_interp-md->eta))) /
(tanh_beta_eta + tanh(md->beta*(1-md->eta)));
}

return (md->beta != 0) ? u_proj : u_interp;
return (md->beta != 0) ? tanh_projection(u_interp, md->beta, md->eta) : u_interp;
}
static void cinterp_tensors(vector3 diag_in_1, cvector3 offdiag_in_1, vector3 diag_in_2,
cvector3 offdiag_in_2, vector3 *diag_out, cvector3 *offdiag_out,
Expand Down Expand Up @@ -2584,8 +2584,8 @@ void material_grids_addgradient_point(double *v, std::complex<double> fields_a,
fashion. Since we aren't checking if each design grid is unique, however,
it's up to the user to only have one unique design grid in this volume.*/
vector3 sz = mg->grid_size;
double *vcur = v, *ucur;
ucur = mg->weights;
double *vcur = v;
double *ucur = mg->weights;
uval = matgrid_val(p, tp, oi, mg);
do {
vector3 pb = to_geom_box_coords(p, &tp->objects[oi]);
Expand All @@ -2599,16 +2599,16 @@ void material_grids_addgradient_point(double *v, std::complex<double> fields_a,
}
// no object tree -- the whole domain is the material grid
if (!tp && is_material_grid(default_material)) {
vector3 pb = to_geom_box_coords(p, &tp->objects[oi]);
map_lattice_coordinates(p.x,p.y,p.z);
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've never used this before, but I assume it's the same as above but assumes a box the size of the simulation cell (the lattice in this case)?

Maybe make sure there's no weird issues with the boundaries (when I read lattice I assume periodic boundaries, which obviously isn't always the case) just make sure the last edge in each dimension is mapped to "1" and not back to "0" (or vice-versa).

(Ideally the two routines use the same mapping logic and the only difference is the default geometry in this case)

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The function map_lattice_coordinates maps elements of the vector3 parameter p into the range [0,1]:

meep/src/meepgeom.cpp

Lines 374 to 381 in 8e2b0a0

void map_lattice_coordinates(double &px, double &py, double &pz) {
px = geometry_lattice.size.x == 0 ? 0
: 0.5 + (px - geometry_center.x) / geometry_lattice.size.x;
py = geometry_lattice.size.y == 0 ? 0
: 0.5 + (py - geometry_center.y) / geometry_lattice.size.y;
pz = geometry_lattice.size.z == 0 ? 0
: 0.5 + (pz - geometry_center.z) / geometry_lattice.size.z;
}

I think you were the one that set up the mapping this way. As long as p is defined to be within geometry_lattice, it should map points along the boundaries correctly.

vector3 sz = mg->grid_size;
double *vcur = v, *ucur;
ucur = mg->weights;
uval = matgrid_val(p, tp, oi, mg);
add_interpolate_weights(pb.x, pb.y, pb.z, vcur, sz.x, sz.y, sz.z, 1,
double *vcur = v;
double *ucur = mg->weights;
uval = material_grid_val(p, mg);
if (mg->beta != 0) uval = tanh_projection(uval, mg->beta, mg->eta);
add_interpolate_weights(p.x, p.y, p.z, vcur, sz.x, sz.y, sz.z, 1,
get_material_gradient(uval, fields_a, fields_f, freq, mg, field_dir) *
scalegrad,
ucur, kind, uval);
tp = geom_tree_search_next(p, tp, &oi);
}
}

Expand Down