diff --git a/python/Makefile.am b/python/Makefile.am index e44fa6216..8df191e1a 100644 --- a/python/Makefile.am +++ b/python/Makefile.am @@ -39,7 +39,9 @@ ADJOINT_TESTS = \ TESTS = \ $(TEST_DIR)/test_3rd_harm_1d.py \ $(TEST_DIR)/test_absorber_1d.py \ + $(TEST_DIR)/test_adjoint_adjacent_grids.py \ $(TEST_DIR)/test_adjoint_chunks.py \ + $(TEST_DIR)/test_adjoint_symmetric_grids.py \ $(TEST_DIR)/test_antenna_radiation.py \ $(TEST_DIR)/test_array_metadata.py \ $(TEST_DIR)/test_bend_flux.py \ diff --git a/python/adjoint/utils.py b/python/adjoint/utils.py index 4ac3abf8f..155af03d4 100644 --- a/python/adjoint/utils.py +++ b/python/adjoint/utils.py @@ -73,6 +73,7 @@ def get_gradient( onp.array(frequencies), sim.geps, finite_difference_step, + getattr(self.design_parameters, "grid_id", -1), ) return onp.squeeze(grad).T diff --git a/python/geom.py b/python/geom.py index a3e54f351..b1207c75b 100755 --- a/python/geom.py +++ b/python/geom.py @@ -7,6 +7,7 @@ from collections import namedtuple from copy import deepcopy import functools +import itertools import math from numbers import Number import operator @@ -585,6 +586,9 @@ def _get_epsmu( return np.squeeze(epsmu) +_material_grid_counter = itertools.count() + + class MaterialGrid: """ This class is used to specify materials on a rectilinear grid. A class object is passed @@ -667,6 +671,12 @@ def __init__( allow you to combine any material grids that overlap in space with no intervening objects. """ self.grid_size = mp.Vector3(*grid_size) + # A stable identity for this grid, carried into the C++ material_data so + # that the adjoint gradient of a design region can tell its own grid + # apart from a neighbouring one. Several material_data objects may be + # created from this single Python object (e.g. one per symmetric copy), + # and they all share this id. + self.grid_id = next(_material_grid_counter) self.medium1 = medium1 self.medium2 = medium2 diff --git a/python/meep.i b/python/meep.i index 9c917969d..c8bdf9c70 100644 --- a/python/meep.i +++ b/python/meep.i @@ -871,7 +871,7 @@ void _get_gradient(PyObject *grad, double scalegrad, meep::dft_fields *fields_a_0, meep::dft_fields *fields_a_1, meep::dft_fields *fields_a_2, meep::dft_fields *fields_f_0, meep::dft_fields *fields_f_1, meep::dft_fields *fields_f_2, meep::grid_volume *grid_volume, PyObject *frequencies, - meep_geom::geom_epsilon *geps, double fd_step) { + meep_geom::geom_epsilon *geps, double fd_step, long owner_grid_id) { // clean the gradient array PyArrayObject *pao_grad = (PyArrayObject *)grad; @@ -896,7 +896,7 @@ void _get_gradient(PyObject *grad, double scalegrad, if (PyArray_DIMS(pao_grad)[0] != nf) meep::abort("Numpy grad array is allocated for %td frequencies; it should be allocated for %td.",PyArray_DIMS(pao_grad)[0],nf); // calculate the gradient - meep_geom::material_grids_addgradient(grad_c,ng,nf,adjoint_fields,forward_fields,frequencies_c,scalegrad,*grid_volume,geps,fd_step); + meep_geom::material_grids_addgradient(grad_c,ng,nf,adjoint_fields,forward_fields,frequencies_c,scalegrad,*grid_volume,geps,owner_grid_id,fd_step); } %} diff --git a/python/tests/test_adjoint_adjacent_grids.py b/python/tests/test_adjoint_adjacent_grids.py new file mode 100644 index 000000000..9acfe036b --- /dev/null +++ b/python/tests/test_adjoint_adjacent_grids.py @@ -0,0 +1,182 @@ +"""Each design region's gradient must come only from its own MaterialGrid. + +Design regions may be adjacent, separated, or overlapping; in every arrangement +the adjoint gradient of one region must agree with a finite difference of the +objective with respect to that region's own design variables. +""" + +import unittest + +import numpy as np +from autograd import numpy as npa + +import meep as mp +import meep.adjoint as mpa + +RESOLUTION = 20 +DPML = 1.0 +CELL = 6.0 +LX, LY = 1.0, 1.0 +FCEN = 1 / 1.55 +SILICON = mp.Medium(index=3.4) +OXIDE = mp.Medium(index=1.44) + + +def _problem(grid_shapes, gap=0.0, overlap=False, grid_type="U_DEFAULT"): + """Design region(s) stacked in y, optionally separated by `gap`. + + With `overlap`, the grids share one volume and only the first is a design + region. + """ + if overlap: + boxes = [(-LY / 2, LY / 2)] * len(grid_shapes) + else: + if len(grid_shapes) == 1: + boxes = [(-LY / 2, LY / 2)] + else: + boxes = [(-LY / 2, -gap / 2), (gap / 2, LY / 2)] + + grids = [ + mp.MaterialGrid( + mp.Vector3(*s), + OXIDE, + SILICON, + # Non-design grids are filled at mid density: an almost-empty cell + # scatters too weakly for the gradient to converge at this decay_by. + weights=np.full(s, 0.5) if overlap else np.zeros(s), + do_averaging=False, + beta=0, + eta=0.5, + grid_type=grid_type, + ) + for s in grid_shapes + ] + vols = [ + mp.Volume(center=mp.Vector3(0, 0.5 * (a + b)), size=mp.Vector3(LX, b - a)) + for a, b in boxes + ] + sim = mp.Simulation( + cell_size=mp.Vector3(CELL, CELL), + resolution=RESOLUTION, + boundary_layers=[mp.PML(DPML)], + geometry=[ + mp.Block(center=v.center, size=v.size, material=g) + for v, g in zip(vols, grids) + ], + sources=[ + mp.Source( + mp.GaussianSource(FCEN, fwidth=0.1), + component=mp.Ez, + center=mp.Vector3(-CELL / 2 + DPML + 0.3, 0), + size=mp.Vector3(0, CELL - 2 * DPML), + ) + ], + default_material=OXIDE, + dimensions=2, + ) + monitor = mp.Volume( + center=mp.Vector3(CELL / 2 - DPML - 0.3, 0), size=mp.Vector3(0, 2.0) + ) + opt = mpa.OptimizationProblem( + simulation=sim, + objective_functions=[lambda ez: npa.sum(npa.abs(ez) ** 2, axis=1)], + objective_arguments=[mpa.FourierFields(sim, monitor, mp.Ez, yee_grid=False)], + design_regions=( + [mpa.DesignRegion(grids[0], volume=vols[0])] + if overlap + else [mpa.DesignRegion(g, volume=v) for g, v in zip(grids, vols)] + ), + frequencies=[FCEN], + decay_by=1e-9, + ) + return opt + + +def _directional_ratio(opt, grid_shapes, seed=0, step=3e-3, n_regions=None): + """adjoint / finite difference along one random direction. + + Directional rather than single-pixel: perturbing one design variable moves + the objective by about as much as the DFT convergence floor. + """ + rng = np.random.default_rng(seed) + if n_regions is None: + n_regions = len(grid_shapes) + weights = [rng.uniform(0.3, 0.7, int(np.prod(s))) for s in grid_shapes[:n_regions]] + _, gradient = opt(rho_vector=[w.copy() for w in weights]) + if n_regions == 1: + gradient = [gradient] + gradient = [np.asarray(g).ravel() for g in gradient] + + def objective(ws): + value, _ = opt(rho_vector=[w.copy() for w in ws], need_gradient=False) + return float(np.ravel(value)[0]) + + direction = [rng.standard_normal(w.size) for w in weights] + norm = np.sqrt(sum(float(np.sum(d**2)) for d in direction)) + direction = [d / norm for d in direction] + + adjoint = sum(float(np.dot(gradient[k], direction[k])) for k in range(len(weights))) + finite_difference = ( + objective([weights[k] + step * direction[k] for k in range(len(weights))]) + - objective([weights[k] - step * direction[k] for k in range(len(weights))]) + ) / (2 * step) + return adjoint / finite_difference + + +class TestAdjacentDesignGrids(unittest.TestCase): + # What the field solve reproduces across the two runs of the central + # difference, not what the gradient assembly contributes. + tol = 2e-3 + + def test_single_design_region(self): + """A single design region covering the whole area.""" + shapes = [(8, 16)] + ratio = _directional_ratio(_problem(shapes), shapes) + self.assertAlmostEqual(ratio, 1.0, delta=self.tol) + + def test_separated_design_regions(self): + """Two design regions separated by a pixel of background.""" + shapes = [(8, 8), (8, 13)] + ratio = _directional_ratio(_problem(shapes, gap=1.0 / RESOLUTION), shapes) + self.assertAlmostEqual(ratio, 1.0, delta=self.tol) + + def test_adjacent_design_regions(self): + """Two design regions sharing a boundary, with equal grid sizes.""" + shapes = [(8, 8), (8, 8)] + ratio = _directional_ratio(_problem(shapes, gap=0.0), shapes) + self.assertAlmostEqual(ratio, 1.0, delta=self.tol) + + def test_adjacent_design_regions_unequal_sizes(self): + """Two design regions sharing a boundary, with unequal grid sizes. + + The larger neighbour's grid_size must not be used to index this + region's smaller gradient array. + """ + shapes = [(8, 8), (8, 13)] + ratio = _directional_ratio(_problem(shapes, gap=0.0), shapes) + self.assertAlmostEqual(ratio, 1.0, delta=self.tol) + + def test_overlapping_grids_are_unaffected(self): + """Two grids overlapping in one volume, combined with U_MEAN.""" + shapes = [(8, 8), (8, 8)] + ratio = _directional_ratio( + _problem(shapes, overlap=True, grid_type="U_MEAN"), shapes, n_regions=1 + ) + self.assertAlmostEqual(ratio, 1.0, delta=self.tol) + + def test_overlapping_grids_unequal_sizes(self): + """An overlapping grid must not supply the owner's interpolation size.""" + # Keep the owner larger than the top grid so the regression is a wrong + # derivative rather than a test that relies on an out-of-bounds access. + shapes = [(8, 13), (8, 8)] + ratio = _directional_ratio( + _problem(shapes, overlap=True, grid_type="U_MEAN"), shapes, n_regions=1 + ) + # Single precision loses accuracy to cancellation in the outer + # objective finite difference; keep the double-precision bound unchanged. + tol = 5e-3 if mp.is_single_precision() else self.tol + self.assertAlmostEqual(ratio, 1.0, delta=tol) + + +if __name__ == "__main__": + unittest.main() diff --git a/python/tests/test_adjoint_symmetric_grids.py b/python/tests/test_adjoint_symmetric_grids.py new file mode 100644 index 000000000..807c06363 --- /dev/null +++ b/python/tests/test_adjoint_symmetric_grids.py @@ -0,0 +1,142 @@ +"""A MaterialGrid shared by several transformed copies of itself. + +Placing one MaterialGrid in several overlapping Blocks, combined with U_MEAN, is +how a design is constrained to a symmetry: a design variable then backs every +copy at once. The adjoint gradient must account for all of them, so it has to +agree with a finite difference taken with respect to that shared variable. +""" + +import unittest + +import numpy as np +from autograd import numpy as npa + +import meep as mp +import meep.adjoint as mpa + +RESOLUTION = 20 +DPML = 1.0 +CELL = 6.0 +LX, LY = 1.0, 1.0 +N = 8 +FCEN = 1 / 1.55 +SILICON = mp.Medium(index=3.4) +OXIDE = mp.Medium(index=1.44) + + +def _problem(transform): + """One design grid, optionally overlapped with a transformed copy of itself. + + `transform` is None, "rot90" or "mirror". The copy reuses the same + MaterialGrid object, which is what makes it the same design variable. + """ + grid = mp.MaterialGrid( + mp.Vector3(N, N), + OXIDE, + SILICON, + weights=np.zeros((N, N)), + do_averaging=False, + beta=0, + eta=0.5, + grid_type="U_MEAN", + ) + volume = mp.Volume(center=mp.Vector3(), size=mp.Vector3(LX, LY)) + blocks = [mp.Block(center=volume.center, size=volume.size, material=grid)] + if transform == "rot90": + blocks.append( + mp.Block( + center=volume.center, + size=volume.size, + material=grid, + e1=mp.Vector3(1).rotate(mp.Vector3(z=1), np.pi / 2), + e2=mp.Vector3(y=1).rotate(mp.Vector3(z=1), np.pi / 2), + ) + ) + elif transform == "mirror": + blocks.append( + mp.Block( + center=volume.center, + size=volume.size, + material=grid, + e1=mp.Vector3(-1, 0, 0), + e2=mp.Vector3(0, 1, 0), + ) + ) + + sim = mp.Simulation( + cell_size=mp.Vector3(CELL, CELL), + resolution=RESOLUTION, + boundary_layers=[mp.PML(DPML)], + geometry=blocks, + sources=[ + mp.Source( + mp.GaussianSource(FCEN, fwidth=0.1), + component=mp.Ez, + center=mp.Vector3(-CELL / 2 + DPML + 0.3, 0), + size=mp.Vector3(0, CELL - 2 * DPML), + ) + ], + default_material=OXIDE, + dimensions=2, + ) + monitor = mp.Volume( + center=mp.Vector3(CELL / 2 - DPML - 0.3, 0), size=mp.Vector3(0, 2.0) + ) + return mpa.OptimizationProblem( + simulation=sim, + objective_functions=[lambda ez: npa.sum(npa.abs(ez) ** 2, axis=1)], + objective_arguments=[mpa.FourierFields(sim, monitor, mp.Ez, yee_grid=False)], + design_regions=[mpa.DesignRegion(grid, volume=volume)], + frequencies=[FCEN], + decay_by=1e-9, + ) + + +def _directional_ratio(opt, seed=0, step=3e-3): + """adjoint / finite difference along one random direction. + + Directional rather than single-pixel: perturbing one design variable moves + the objective by about as much as the DFT convergence floor. + """ + rng = np.random.default_rng(seed) + weights = rng.uniform(0.3, 0.7, N * N) + _, gradient = opt(rho_vector=[weights.copy()]) + gradient = np.asarray(gradient).ravel() + + def objective(w): + value, _ = opt(rho_vector=[w.copy()], need_gradient=False) + return float(np.ravel(value)[0]) + + direction = rng.standard_normal(weights.size) + direction /= np.linalg.norm(direction) + adjoint = float(np.dot(gradient, direction)) + finite_difference = ( + objective(weights + step * direction) - objective(weights - step * direction) + ) / (2 * step) + return adjoint / finite_difference + + +class TestSymmetricDesignGrids(unittest.TestCase): + # What the field solve reproduces across the two runs of the central + # difference, not what the gradient assembly contributes. + tol = 2e-3 + + def test_no_transform(self): + """A single block, with no overlapping copy.""" + self.assertAlmostEqual(_directional_ratio(_problem(None)), 1.0, delta=self.tol) + + def test_rotated_copy(self): + """Overlapped with a 90-degree rotation of itself.""" + self.assertAlmostEqual( + _directional_ratio(_problem("rot90")), 1.0, delta=self.tol + ) + + def test_mirrored_copy(self): + """Overlapped with a reflection of itself.""" + self.assertAlmostEqual( + _directional_ratio(_problem("mirror")), 1.0, delta=self.tol + ) + + +if __name__ == "__main__": + unittest.main() diff --git a/python/typemap_utils.cpp b/python/typemap_utils.cpp index e56dff4fc..622151262 100644 --- a/python/typemap_utils.cpp +++ b/python/typemap_utils.cpp @@ -495,6 +495,16 @@ static int pymaterial_grid_to_material_grid(PyObject *po, material_data *md) { // initialize grid size if (!get_attr_v3(po, &md->grid_size, "grid_size")) { return 0; } + // carry the Python-level identity across, so the adjoint gradient of a design + // region can tell its own grid apart from a neighbouring one + { + PyObject *py_grid_id = PyObject_GetAttrString(po, "grid_id"); + if (!py_grid_id) return 0; + md->grid_id = PyLong_AsLong(py_grid_id); + Py_DECREF(py_grid_id); + if (md->grid_id == -1 && PyErr_Occurred()) return 0; + } + // initialize user specified materials PyObject *po_medium1 = PyObject_GetAttrString(po, "medium1"); if (!po_medium1 || !pymedium_to_medium(po_medium1, &md->medium_1)) { diff --git a/src/material_data.cpp b/src/material_data.cpp index 2b2f1cfd9..f675d4602 100644 --- a/src/material_data.cpp +++ b/src/material_data.cpp @@ -46,8 +46,8 @@ void medium_struct::check_offdiag_im_zero_or_abort() const { material_data::material_data() : which_subclass(MEDIUM), medium(), user_func(NULL), user_data(NULL), do_averaging(false), - epsilon_data(NULL), epsilon_dims{}, grid_size{}, weights(NULL), medium_1(), medium_2(), - material_grid_kinds{U_DEFAULT} {} + epsilon_data(NULL), epsilon_dims{}, grid_size{}, weights(NULL), grid_id(-1), medium_1(), + medium_2(), material_grid_kinds{U_DEFAULT} {} void material_data::copy_from(const material_data &from) { which_subclass = from.which_subclass; @@ -72,6 +72,8 @@ void material_data::copy_from(const material_data &from) { std::copy_n(from.weights, N, weights); } + grid_id = from.grid_id; + medium_1 = from.medium_1; medium_2 = from.medium_2; beta = from.beta; diff --git a/src/material_data.hpp b/src/material_data.hpp index 28f803b47..417d5fff9 100644 --- a/src/material_data.hpp +++ b/src/material_data.hpp @@ -131,6 +131,7 @@ struct material_data { // these fields used only if which_subclass==MATERIAL_GRID vector3 grid_size; double *weights; + long grid_id; medium_struct medium_1; medium_struct medium_2; double beta; diff --git a/src/meepgeom.cpp b/src/meepgeom.cpp index 3bd92f875..dce99195e 100644 --- a/src/meepgeom.cpp +++ b/src/meepgeom.cpp @@ -15,6 +15,7 @@ */ #include +#include #include #include "meepgeom.hpp" #include "meep_internals.hpp" @@ -2640,7 +2641,8 @@ get_material_gradient(const meep::vec &r, // current point geom_epsilon *geps, // material meep::grid_volume &gv, // simulation grid volume double du, // step size - double *u, // matgrid + const std::vector &us, // matgrid weights, one per + // copy of this design grid int idx // matgrid index ) { /*Compute the Aᵤx product from the -λᵀAᵤx calculation. @@ -2686,23 +2688,36 @@ get_material_gradient(const meep::vec &r, // current point v.set_direction_max(d, r.in_direction(d) + 0.5 * gv.inva * sd); } double row_1[3], row_2[3]; - double orig = u[idx]; - u[idx] -= du; + /* One design variable may back several copies of the same grid (the way a + symmetry is imposed), so it has to be perturbed in every copy at once -- + perturbing one copy answers a different question. */ + std::vector orig(us.size()); + for (size_t k = 0; k < us.size(); ++k) { + orig[k] = us[k][idx]; + us[k][idx] -= du; + } geps->eff_chi1inv_row(adjoint_c, row_1, v, geps->tol, geps->maxeval); - u[idx] += 2 * du; + for (size_t k = 0; k < us.size(); ++k) + us[k][idx] = orig[k] + du; geps->eff_chi1inv_row(adjoint_c, row_2, v, geps->tol, geps->maxeval); - u[idx] = orig; + for (size_t k = 0; k < us.size(); ++k) + us[k][idx] = orig[k]; return fields_f * (row_1[dir_idx] - row_2[dir_idx]) / (2 * du); } // materials have some dispersion else { - double orig = u[idx]; + std::vector orig(us.size()); std::complex row_1[3], row_2[3], dA_du[3]; - u[idx] -= du; + for (size_t k = 0; k < us.size(); ++k) { + orig[k] = us[k][idx]; + us[k][idx] -= du; + } eff_chi1inv_row_disp(adjoint_c, row_1, r, freq, geps); - u[idx] += 2 * du; + for (size_t k = 0; k < us.size(); ++k) + us[k][idx] = orig[k] + du; eff_chi1inv_row_disp(adjoint_c, row_2, r, freq, geps); - u[idx] = orig; + for (size_t k = 0; k < us.size(); ++k) + us[k][idx] = orig[k]; return fields_f * (row_1[dir_idx] - row_2[dir_idx]) / (2 * du) * cond_cmp(forward_c, r, freq, geps); @@ -2716,58 +2731,66 @@ With the addition of subpixel smoothing, however, the vJp became much more complicated and it is easier to calculate the entire gradient using finite differences (at the cost of slightly less accurate gradients due to floating-point roundoff errors). */ -void add_interpolate_weights(double rx, double ry, double rz, double *data, int nx, int ny, int nz, - int stride, double scaleby, double *udata, int ukind, double uval, - meep::vec r, geom_epsilon *geps, meep::component adjoint_c, - meep::component forward_c, std::complex fwd, - std::complex adj, double freq, meep::grid_volume &gv, - double du) { - int x1, y1, z1, x2, y2, z2; - double dx, dy, dz, u; +void add_interpolate_weights(const std::vector &udatas, const std::vector &pbs, + double *data, int nx, int ny, int nz, int stride, double scaleby, + int ukind, double uval, meep::vec r, geom_epsilon *geps, + meep::component adjoint_c, meep::component forward_c, + std::complex fwd, std::complex adj, double freq, + meep::grid_volume &gv, double du) { + /* `udatas`/`pbs` are the copies of *one* design grid that cover this point, + with the box coordinates of each. A design variable is shared across all of + them, so it is differentiated once, with every copy perturbed together, and + the union of the copies' interpolation stencils is the set of variables this + point touches. */ - meep::map_coordinates(rx, ry, rz, nx, ny, nz, x1, y1, z1, x2, y2, z2, dx, dy, dz); - int x_list[2] = {x1, x2}, y_list[2] = {y1, y2}, z_list[2] = {z1, z2}; - int lx = (x1 == x2) ? 1 : 2; - int ly = (y1 == y2) ? 1 : 2; - int lz = (z1 == z2) ? 1 : 2; - -/* define a macro to give us data(x,y,z) on the grid, -in row-major order (the order used by HDF5): */ #define IDX(x, y, z) (((x)*ny + (y)) * nz + (z)) * stride -#define D(x, y, z) (data[IDX(x, y, z)]) -#define U(x, y, z) (udata[IDX(x, y, z)]) - u = (((U(x1, y1, z1) * (1.0 - dx) + U(x2, y1, z1) * dx) * (1.0 - dy) + - (U(x1, y2, z1) * (1.0 - dx) + U(x2, y2, z1) * dx) * dy) * - (1.0 - dz) + - ((U(x1, y1, z2) * (1.0 - dx) + U(x2, y1, z2) * dx) * (1.0 - dy) + - (U(x1, y2, z2) * (1.0 - dx) + U(x2, y2, z2) * dx) * dy) * - dz); - - if (ukind == material_data::U_MIN && u != uval) return; // TODO look into this - if (ukind == material_data::U_PROD) scaleby *= uval / u; - - for (int xi = 0; xi < lx; xi++) { - for (int yi = 0; yi < ly; yi++) { - for (int zi = 0; zi < lz; zi++) { - int x = x_list[xi], y = y_list[yi], z = z_list[zi]; - int u_idx = IDX(x, y, z); - std::complex prod = adj * get_material_gradient(r, adjoint_c, forward_c, fwd, freq, - geps, gv, du, udata, u_idx); - D(x, y, z) += prod.real() * scaleby; - } + std::set touched; + for (size_t k = 0; k < pbs.size(); ++k) { + int x1, y1, z1, x2, y2, z2; + double dx, dy, dz; + meep::map_coordinates(pbs[k].x, pbs[k].y, pbs[k].z, nx, ny, nz, x1, y1, z1, x2, y2, z2, dx, dy, + dz); + int x_list[2] = {x1, x2}, y_list[2] = {y1, y2}, z_list[2] = {z1, z2}; + int lx = (x1 == x2) ? 1 : 2; + int ly = (y1 == y2) ? 1 : 2; + int lz = (z1 == z2) ? 1 : 2; + + if (k == 0 && (ukind == material_data::U_MIN || ukind == material_data::U_PROD)) { + /* U_MIN and U_PROD carry adjustments defined in terms of the interpolated + value at this point. Overlapping grids of these kinds are rejected by + the caller, so there is exactly one copy and this is unchanged. */ + const double *U = udatas[k]; + double u = (((U[IDX(x1, y1, z1)] * (1.0 - dx) + U[IDX(x2, y1, z1)] * dx) * (1.0 - dy) + + (U[IDX(x1, y2, z1)] * (1.0 - dx) + U[IDX(x2, y2, z1)] * dx) * dy) * + (1.0 - dz) + + ((U[IDX(x1, y1, z2)] * (1.0 - dx) + U[IDX(x2, y1, z2)] * dx) * (1.0 - dy) + + (U[IDX(x1, y2, z2)] * (1.0 - dx) + U[IDX(x2, y2, z2)] * dx) * dy) * + dz); + if (ukind == material_data::U_MIN && u != uval) return; // TODO look into this + if (ukind == material_data::U_PROD) scaleby *= uval / u; } + + for (int xi = 0; xi < lx; xi++) + for (int yi = 0; yi < ly; yi++) + for (int zi = 0; zi < lz; zi++) + touched.insert(IDX(x_list[xi], y_list[yi], z_list[zi])); + } + + for (std::set::const_iterator it = touched.begin(); it != touched.end(); ++it) { + std::complex prod = + adj * get_material_gradient(r, adjoint_c, forward_c, fwd, freq, geps, gv, du, udatas, *it); + data[*it] += prod.real() * scaleby; } #undef IDX -#undef D -#undef U } void material_grids_addgradient_point(double *v, vector3 p, double scalegrad, geom_epsilon *geps, meep::component adjoint_c, meep::component forward_c, std::complex fwd, std::complex adj, - double freq, meep::grid_volume &gv, double tol) { + double freq, meep::grid_volume &gv, double tol, + long owner_grid_id) { geom_box_tree tp; int oi, ois; material_data *mg, *mg_sum; @@ -2794,7 +2817,10 @@ void material_grids_addgradient_point(double *v, vector3 p, double scalegrad, ge ++matgrid_val_count; if (tp_sum) mg_sum = (material_data *)tp_sum->objects[ois].o->material; } while (tp_sum && is_material_grid(mg_sum)); - scalegrad /= matgrid_val_count; + /* No scalegrad /= matgrid_val_count here: the finite difference in + get_material_gradient() runs through matgrid_val(), so the 1/N of the + U_MEAN average is already contained in the result. */ + (void)matgrid_val_count; } else if ((tp) && ((mg->material_grid_kinds == material_data::U_MIN) || (mg->material_grid_kinds == material_data::U_PROD))) { @@ -2817,26 +2843,47 @@ void material_grids_addgradient_point(double *v, vector3 p, double scalegrad, ge 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; - double *ucur = mg->weights; uval = tanh_projection(matgrid_val(p, tp, oi, mg), mg->beta, mg->eta); + + /* Collect the copies of this design region's own grid that cover the point. + + Only the owner's copies: v[] is sized and indexed by the owner's + parameterization, so a contribution from any other grid is mis-attributed, + and an out-of-bounds write when that grid is larger. That happens whenever + two grids abut, because the design region's DFT monitor is snapped outward + to enclosing grid nodes and so covers nodes shared with the neighbour. The + neighbouring design region visits them on its own pass, so nothing is + lost. All copies of one grid share a grid_id, so a grid overlapped with + transformed copies of itself -- the documented way to impose a symmetry -- + still contributes through every copy. */ + std::vector udatas; + std::vector pbs; + geom_box_tree tpc = tp; + int oic = oi; do { - vector3 pb = to_geom_box_coords(p, &tp->objects[oi]); - add_interpolate_weights(pb.x, pb.y, pb.z, vcur, sz.x, sz.y, sz.z, 1, scalegrad, ucur, kind, - uval, vector3_to_vec(p), geps, adjoint_c, forward_c, fwd, adj, freq, - gv, tol); + material_data *mg_cur = (material_data *)tpc->objects[oic].o->material; + if (owner_grid_id < 0 || mg_cur->grid_id == owner_grid_id) { + if (udatas.empty()) sz = mg_cur->grid_size; + udatas.push_back(mg_cur->weights); + pbs.push_back(to_geom_box_coords(p, &tpc->objects[oic])); + } if (kind == material_data::U_DEFAULT) break; - tp = geom_tree_search_next(p, tp, &oi); - } while (tp && is_material_grid((material_data *)tp->objects[oi].o->material)); + tpc = geom_tree_search_next(p, tpc, &oic); + } while (tpc && is_material_grid((material_data *)tpc->objects[oic].o->material)); + + if (!udatas.empty()) + add_interpolate_weights(udatas, pbs, v, sz.x, sz.y, sz.z, 1, scalegrad, kind, uval, + vector3_to_vec(p), geps, adjoint_c, forward_c, fwd, adj, freq, gv, + tol); } // no object tree -- the whole domain is the material grid if (!tp && is_material_grid(default_material)) { map_lattice_coordinates(p.x, p.y, p.z); vector3 sz = mg->grid_size; - double *vcur = v; - double *ucur = mg->weights; uval = tanh_projection(material_grid_val(p, mg), mg->beta, mg->eta); - add_interpolate_weights(p.x, p.y, p.z, vcur, sz.x, sz.y, sz.z, 1, scalegrad, ucur, kind, uval, + std::vector udatas(1, mg->weights); + std::vector pbs(1, p); + add_interpolate_weights(udatas, pbs, v, sz.x, sz.y, sz.z, 1, scalegrad, kind, uval, vector3_to_vec(p), geps, adjoint_c, forward_c, fwd, adj, freq, gv, tol); } } @@ -2899,7 +2946,7 @@ void material_grids_addgradient(double *v, size_t ng, size_t nf, std::vector fields_a, std::vector fields_f, double *frequencies, double scalegrad, meep::grid_volume &gv, geom_epsilon *geps, - double du) { + long owner_grid_id, double du) { /* ------------------------------------------------------------ */ // initialize local gradient array /* ------------------------------------------------------------ */ @@ -2997,7 +3044,7 @@ void material_grids_addgradient(double *v, size_t ng, size_t nf, : 1; // the pi is already factored in near2far.cpp material_grids_addgradient_point(v_local + ng * f_i, vec_to_vector3(p), scalegrad * cyl_scale, geps, adjoint_c, forward_c, - fwd, adj, frequencies[f_i], gv, du); + fwd, adj, frequencies[f_i], gv, du, owner_grid_id); /* more complicated case requires interpolation/restriction */ } else if ((md->do_averaging) || /* account for subpixel smoothing */ @@ -3039,7 +3086,7 @@ void material_grids_addgradient(double *v, size_t ng, size_t nf, material_grids_addgradient_point(v_local + ng * f_i, vec_to_vector3(eps1), scalegrad * cyl_scale, geps, adjoint_c, forward_c, fwd_avg, std::complex(0.5, 0) * adj, - frequencies[f_i], gv, du); + frequencies[f_i], gv, du, owner_grid_id); } } /********* compute λᵀbᵤ ***************/ diff --git a/src/meepgeom.hpp b/src/meepgeom.hpp index b2a0d3e5b..004954209 100644 --- a/src/meepgeom.hpp +++ b/src/meepgeom.hpp @@ -298,7 +298,7 @@ void material_grids_addgradient(double *v, size_t ng, size_t nf, std::vector fields_a, std::vector fields_f, double *frequencies, double scalegrad, meep::grid_volume &gv, geom_epsilon *geps, - double du = 1e-6); + long owner_grid_id, double du = 1e-6); /***************************************************************/ /* routines in GDSIIgeom.cc ************************************/