diff --git a/NEWS.md b/NEWS.md index 0bcde63bd..21d36a170 100644 --- a/NEWS.md +++ b/NEWS.md @@ -2,6 +2,17 @@ ## Meep 1.35.0 (in progress) +* Adjoint solver: geometric objects can be differentiated with respect to their + `center` and `size`, so an optimizer can be asked whether a reflector should + move or a spacer should lengthen. An object opts in with + `differentiable=['center', 'size']` and its gradient appears under its `name`. + It costs no extra simulation. Only pixels the object's boundary passes through + contribute, which is the discrete form of a shape derivative being a surface + integral. Requires subpixel smoothing, and refuses to run without it. Object + faces should be kept a quarter pixel clear of pixel edges, since Yee + components sit half a pixel apart and a face on some component's pixel edge + has only a one-sided derivative; a warning is issued when this is detected. + * Adjoint solver: sources can now be differentiated alongside the design regions. A source opts in by naming parameters, as in `differentiable=['beam_w0', 'beam_x0']`, and its gradient appears in the diff --git a/doc/docs/Python_Tutorials/Adjoint_Solver.md b/doc/docs/Python_Tutorials/Adjoint_Solver.md index b5245afe1..86c47f87a 100644 --- a/doc/docs/Python_Tutorials/Adjoint_Solver.md +++ b/doc/docs/Python_Tutorials/Adjoint_Solver.md @@ -492,6 +492,72 @@ the run was doubled. Finally, a source inside or near a PML is rejected: its adjoint field is absorbed, so the gradient would come back finite, smooth, and wrong. +Differentiating With Respect To Geometry +---------------------------------------- + +The design region is not the only thing that can move. An ordinary geometric +object can be differentiated too — where a reflector sits, how long a spacer is +— by naming the parameters on the object itself: + +```py +reflector = mp.Block( + center=mp.Vector3(0, 2.0), + size=mp.Vector3(4, 0.3), + material=silicon, + differentiable=["center", "size"], + name="reflector", +) + +value, grad = opt([rho]) +grad["design"] # as before +grad["reflector"]["center"] # (3,), or (nfreq, 3) for several frequencies +grad["reflector"]["size"] +``` + +This is the same pattern as the design weights and the Gaussian beam +parameters: finite-difference a cheap analytic map and contract it against the +adjoint field, at no cost in extra timestepping. Subpixel smoothing makes the +permittivity depend on the geometry only through each pixel's filling fraction +and the interface normal, so + +$$\frac{\partial \chi^{-1}}{\partial p} = \frac{\partial \chi^{-1}}{\partial f}\,\frac{\partial f}{\partial p}$$ + +with $\partial f/\partial p$ analytic — a block is an intersection of slabs, so +the pixel overlap factorizes — and only pixels the boundary passes through +contributing anything. That is the discrete form of a shape derivative being a +surface integral. + +`'center'` and `'size'` are accepted here and rejected on a *source*. That is +not an inconsistency: a source's cotangent is gathered over a fixed set of grid +points, so moving it changes which points it occupies rather than the +amplitudes applied to them, whereas the permittivity is a function of position +and moving an object is exactly what a derivative with respect to position +means. + +### Two things to get right + +**Subpixel smoothing must be on.** It is on by default (`eps_averaging=True`), +and the gradient refuses to run without it rather than returning a number. +Without smoothing the permittivity is a step function of position — nothing +changes until a boundary crosses a pixel edge, then it changes by the full +material contrast — and a finite difference of that is not a derivative. + +**Keep object faces off the pixel edges.** Pixel centres lie at integer +multiples of the pixel, so pixel edges lie at half-integers, and Yee components +sit half a pixel apart from one another. A face at a pixel centre for one field +component therefore lies exactly on a pixel edge for another, where no pixel +straddles it and the derivative is one-sided. A **quarter-pixel** offset clears +both, so every component straddles every face: + +```py +offset = 0.25 / resolution +block = mp.Block(center=mp.Vector3(0, y0 + offset), ...) +``` + +Round geometry on a round grid lands on edges constantly — a block 1.0 wide at +resolution 20 has faces exactly 10 pixels from its centre — so this is worth +doing deliberately. A warning is issued when a face is detected on an edge. + Broadband Waveguide Mode Converter with Minimum Feature Size ------------------------------------------------------------ diff --git a/python/Makefile.am b/python/Makefile.am index 985c3df58..efe5e4ff3 100644 --- a/python/Makefile.am +++ b/python/Makefile.am @@ -37,6 +37,7 @@ ADJOINT_TESTS = \ $(TEST_DIR)/test_adjoint_symmetry.py \ $(TEST_DIR)/test_adjoint_protocol.py \ $(TEST_DIR)/test_angular_spectrum.py \ + $(TEST_DIR)/test_geometry_gradient.py \ $(TEST_DIR)/test_source_gradient.py \ $(TEST_DIR)/test_adjoint_jax.py @@ -101,6 +102,7 @@ TESTS = \ $(TEST_DIR)/test_simulation.py \ $(TEST_DIR)/test_special_kz.py \ $(TEST_DIR)/test_source.py \ + $(TEST_DIR)/test_geometry_gradient.py \ $(TEST_DIR)/test_source_gradient.py \ $(TEST_DIR)/test_stop_when_flux_decayed.py \ $(TEST_DIR)/test_subpixel_3d.py \ @@ -256,6 +258,7 @@ adjoint_PYTHON = $(srcdir)/adjoint/__init__.py \ $(srcdir)/adjoint/filter_source.py \ $(srcdir)/adjoint/connectivity.py \ $(srcdir)/adjoint/unfilter_design.py \ + $(srcdir)/adjoint/geometry_gradient.py \ $(srcdir)/adjoint/source_gradient.py \ $(srcdir)/adjoint/wrapper.py \ $(srcdir)/adjoint/utils.py diff --git a/python/adjoint/__init__.py b/python/adjoint/__init__.py index 16c1b3fdd..d71c614cc 100644 --- a/python/adjoint/__init__.py +++ b/python/adjoint/__init__.py @@ -21,6 +21,7 @@ from .unfilter_design import * from . import source_gradient +from . import geometry_gradient # JAX is an optional dependency; everything that needs it lives in `wrapper`. # Importing it also registers JAX as a way to differentiate objective functions, diff --git a/python/adjoint/geometry_gradient.py b/python/adjoint/geometry_gradient.py new file mode 100644 index 000000000..d35490a25 --- /dev/null +++ b/python/adjoint/geometry_gradient.py @@ -0,0 +1,295 @@ +"""Adjoint gradients with respect to a geometric object's centre and size. + +The design gradient finite-differences `eff_chi1inv_row` -- the subpixel +smoothed inverse permittivity over a voxel -- with respect to a design weight. +That function is a function of the *geometry*, so perturbing an object's centre +or size instead gives the shape derivative from the same machinery, at no cost +in additional timestepping. It is the same pattern as the design weights and as +the Gaussian beam parameters: finite-difference a cheap analytic map, contract +against the adjoint field. + +Two things differ from the design gradient. + +A design weight is local, so it can be perturbed inside the point loop. A centre +or a size is global, so the perturbation hoists out: move the geometry once, +visit every point, restore. Twelve passes for six parameters, not twelve +perturbations per point. + +And the step has units. `utils.FD_DEFAULT` is a dimensionless perturbation of a +weight; here it is a length, and the scale that makes sense is a fraction of a +pixel. `DEFAULT_STEP_PIXELS` sets it relative to the grid rather than absolutely. + +Accuracy: a shape derivative is one order worse than the fields +------------------------------------------------------------------ +Subpixel smoothing makes the permittivity -- and so the objective -- a +second-order accurate, *continuous* function of an object's position. Measured +on a block of index 2.5 at resolution 20, turning smoothing off makes the +objective piecewise constant in position, swinging 77% across a single pixel; +turning it on makes it continuous, swinging 14%. + +But second-order accuracy in the value is only first-order accuracy in the +derivative. Writing the discretization error as a function of the sub-pixel +phase, + + J(p) = J_exact(p) + Delta^2 E(p / Delta) + +with E an O(1) oscillatory function -- which is why the residual has a period of +exactly one pixel -- differentiating gives + + dJ/dp = dJ_exact/dp + Delta^2 (1/Delta) E'(p / Delta) + = dJ_exact/dp + O(Delta) + +Measured: the per-pixel swing falls 9.9% -> 4.0% -> 1.9% at resolutions +20 -> 40 -> 80, halving with each doubling rather than quartering. This is +inherent to differentiating a smoothed staircase and cannot be recovered by +improving the adjoint, because the error is in the function being +differentiated, not in how it is differentiated. + +Two consequences for testing. The adjoint should be checked against a finite +difference of the *discrete* objective, since both differentiate the same +smoothed function and should agree closely; that is a statement about the +implementation. Agreement with the physically intended derivative is limited to +O(Delta) and is a statement about the discretization. +""" + +import warnings +from typing import List, Optional + +import numpy as np + +import meep as mp + +# Perturbation for the geometric finite difference, as a fraction of a pixel. +# An absolute step is the wrong idea: subpixel smoothing varies over a pixel, so +# the step has to be defined against the grid. +# +# Measured on a block whose faces sit mid-voxel, against a finite difference of +# the objective: 0.2 px is 10% out, 0.05 px is 6.5%, 0.02 px is 3.4%, 0.005 px +# is 1.3% and 0.001 px is 0.7%, still improving. The error is truncation, not +# roundoff, over that whole range, so the step wants to be small. +DEFAULT_STEP_PIXELS = 0.002 + +# How far past an object the gradient monitor reaches, in pixels. The support of +# d(epsilon)/d(parameter) is the shell of voxels the boundary sweeps through, so +# it extends outside the object itself. +MONITOR_PAD_PIXELS = 2.0 + +# Indices `geometry_addgradient` uses, matching the geom_param enumeration in +# meepgeom.cpp. +PARAMETER_INDEX = { + "center": (0, 1, 2), + "size": (3, 4, 5), +} + +DIFFERENTIABLE_PARAMS = tuple(PARAMETER_INDEX) + + +def differentiable_objects(sim: mp.Simulation) -> List: + """The geometric objects flagged for differentiation, with their indices.""" + return [ + (index, obj) + for index, obj in enumerate(sim.geometry) + if getattr(obj, "differentiable", ()) + ] + + +def object_key(obj, index: int): + """The key an object's gradient appears under: its name, or its position.""" + name = getattr(obj, "name", None) + return name if name is not None else index + + +def _check_supported(obj) -> None: + """Refuse the cases the diagonal-only contraction would get wrong.""" + if not isinstance(obj, mp.Block): + raise NotImplementedError( + f"Geometry gradients are implemented for mp.Block, not " + f"{type(obj).__name__}. 'center' is meaningful for any object, but " + "'size' is not, and the contraction currently assumes axis-aligned " + "faces." + ) + axes = (obj.e1, obj.e2, obj.e3) + expected = (mp.Vector3(1, 0, 0), mp.Vector3(0, 1, 0), mp.Vector3(0, 0, 1)) + for axis, unit in zip(axes, expected): + if abs(axis.x - unit.x) + abs(axis.y - unit.y) + abs(axis.z - unit.z) > 1e-12: + raise NotImplementedError( + "Geometry gradients assume an axis-aligned block, so that " + "'size' means what the three components suggest. The " + "off-diagonal terms subpixel smoothing produces at a tilted " + "boundary *are* contracted; it is the parameterization that is " + "restricted, not the physics." + ) + + +def check_faces_off_voxel_edges(sim: mp.Simulation, obj) -> None: + """Warn when a face sits exactly on a voxel edge, where no derivative exists. + + A voxel's filling fraction is piecewise linear in the position of the + boundary crossing it, with a kink each time the boundary reaches a voxel + edge. The smoothed permittivity -- and so the objective -- is therefore C0 + but not C1 in an object's position, and *on* an edge the left and right + derivatives differ. + + A central difference straddles that kink and returns a weighted mixture of + the two, with weights that depend on the step. Measured on such a face, the + reported gradient swept monotonically from -0.19 to +0.12 as the step went + from 0.2 px to 0.001 px, passing through the true value without settling. + Moving the same face half a pixel made the sweep converge to 0.7%. + + This is easy to hit by accident, because round sizes at round resolutions + land on edges: a 0.8-wide block at resolution 20 has faces exactly 8 pixels + from its centre. + """ + dx = 1.0 / sim.resolution + size = (obj.size.x, obj.size.y, obj.size.z) + center = (obj.center.x, obj.center.y, obj.center.z) + for axis, letter in enumerate("xyz"): + if size[axis] == 0: + continue + for face in (center[axis] - size[axis] / 2, center[axis] + size[axis] / 2): + # Pixel centres sit at integer multiples of dx, so pixel *edges* + # are at half-integers. A face at a pixel centre is straddled by + # that pixel and smooths normally; a face at a pixel edge is + # straddled by nothing, every neighbouring pixel is full or empty, + # and the derivative there is one-sided. + offset = abs((face / dx) - round(face / dx)) + if abs(offset - 0.5) < 1e-6: + warnings.warn( + f"The {letter} face of " + f"{getattr(obj, 'name', None) or 'this object'} at " + f"{face:.6g} lies on a pixel edge, where the smoothed " + "permittivity has a kink and the derivative with respect " + "to position is one-sided. The reported gradient will " + "depend on the finite-difference step rather than " + f"converging. Shift it by about {dx / 2:.6g} along " + + letter + + " so the face lands on a pixel centre, or change the " + "resolution.", + RuntimeWarning, + stacklevel=3, + ) + + +def check_smoothing(sim: mp.Simulation) -> None: + """Refuse to differentiate a staircase. + + Without subpixel smoothing the effective permittivity is a step function of + an object's position: shifting a boundary changes nothing until it crosses a + voxel edge, then changes by the full material contrast. A finite difference + of that is zero almost everywhere and a spike in a few voxels, and the sum + is not a derivative. It would look like a plausible number rather than an + error, so refuse instead. + """ + if not getattr(sim, "eps_averaging", True): + raise ValueError( + "Geometry gradients need subpixel smoothing, but this simulation " + "has eps_averaging=False. Without it the permittivity is a step " + "function of an object's position and the derivative does not " + "exist at the grid scale." + ) + + +def install_geometry_monitors( + sim: mp.Simulation, + objects, + frequencies, + decimation_factor: int = 0, +) -> List[List[mp.DftFields]]: + """DFT monitors over each differentiable object's support. + + Unlike a design region, an ordinary object has no monitor of its own, so the + fields it needs have to be recorded in both the forward and adjoint runs. + """ + from . import utils + + monitors = [] + for _, obj in objects: + # Pad beyond the object. Moving a boundary changes the smoothed + # permittivity in every voxel the boundary passes through, and that + # includes voxels whose centres lie *outside* the original extent. A + # monitor covering only the object misses them, which costs little + # where the field is continuous across the boundary and a great deal + # where it is not -- the normal E field jumps by the index contrast. + pad = MONITOR_PAD_PIXELS / sim.resolution + + # `mp.inf` is how one writes "spans the cell", and is the natural way to + # describe a layer -- but it is 1e20, not a flag, and + # `_fit_volume_to_simulation` passes it straight through. A DFT volume + # of that extent fails inside meep with "impossible(?) looping + # boundaries", so clamp to the cell before fitting. + cell = sim.cell_size + padded = mp.Vector3( + *( + min(s + 2 * pad, c) if s else 0.0 + for s, c in zip( + (obj.size.x, obj.size.y, obj.size.z), (cell.x, cell.y, cell.z) + ) + ) + ) + volume = sim._fit_volume_to_simulation( + mp.Volume(center=obj.center, size=padded) + ) + monitors.append( + [ + sim.add_dft_fields( + [component], + frequencies, + where=volume, + yee_grid=True, + decimation_factor=decimation_factor, + persist=True, + ) + for component in utils._compute_components(sim) + ] + ) + return monitors + + +def gradient( + sim: mp.Simulation, + obj, + object_index: int, + forward_fields: List[mp.DftFields], + adjoint_fields: List[mp.DftFields], + frequencies, + step: Optional[float] = None, +) -> dict: + """dJ/d(parameter) for one object, keyed by the names it declared.""" + _check_supported(obj) + check_smoothing(sim) + check_faces_off_voxel_edges(sim, obj) + + names = [n for n in obj.differentiable] + indices = [] + for name in names: + indices.extend(PARAMETER_INDEX[name]) + indices = np.asarray(indices, dtype=np.intc) + + frequencies = np.asarray(frequencies, dtype=np.float64) + out = np.zeros((frequencies.size, indices.size), dtype=np.float64) + + if step is None: + step = DEFAULT_STEP_PIXELS / sim.resolution + + mp._get_geometry_gradient( + out, + 1.0, + adjoint_fields[0].swigobj, + adjoint_fields[1].swigobj, + adjoint_fields[2].swigobj, + forward_fields[0].swigobj, + forward_fields[1].swigobj, + forward_fields[2].swigobj, + sim.gv, + frequencies, + sim.geps, + object_index, + indices, + step, + ) + + grads = {} + for position, name in enumerate(names): + columns = slice(3 * position, 3 * position + 3) + grads[name] = np.squeeze(out[:, columns]) + return grads diff --git a/python/adjoint/optimization_problem.py b/python/adjoint/optimization_problem.py index dc6a38097..36b387e25 100644 --- a/python/adjoint/optimization_problem.py +++ b/python/adjoint/optimization_problem.py @@ -6,6 +6,7 @@ from . import LDOS, DesignRegion, utils, ObjectiveQuantity from . import source_gradient +from . import geometry_gradient class OptimizationProblem: @@ -154,6 +155,12 @@ def __init__( self.differentiable_sources = source_gradient.differentiable_sources(self.sim) self.source_gradient = {} + # Geometric objects flagged with `differentiable=[...]`: the shape + # derivative of the structure itself, as opposed to the density inside + # a design region. + self.differentiable_objects = geometry_gradient.differentiable_objects(self.sim) + self.geometry_gradient = {} + # The optimizer has three allowable states : "INIT", "FWD", and "ADJ". # INIT - The optimizer is initialized and ready to run a forward simulation # FWD - The optimizer has already run a forward simulation @@ -227,10 +234,14 @@ def __call__( f"Incorrect solver state detected: {self.current_state}" ) - if self.differentiable_sources: - # Only change the return shape when the user asked for source - # gradients; without a flagged source this is exactly as before. - return self.f0, {"design": self.gradient, **self.source_gradient} + if self.differentiable_sources or self.differentiable_objects: + # Only change the return shape when something extra was flagged; + # otherwise this is exactly as before. + return self.f0, { + "design": self.gradient, + **self.source_gradient, + **self.geometry_gradient, + } return self.f0, self.gradient @@ -270,6 +281,14 @@ def prepare_forward_run(self): self.forward_design_region_monitors = utils.install_design_region_monitors( self.sim, self.design_regions, self.frequencies, self.decimation_factor ) + # an ordinary object has no monitor of its own, so its fields have to be + # recorded in both runs + self.forward_geometry_monitors = geometry_gradient.install_geometry_monitors( + self.sim, + self.differentiable_objects, + self.frequencies, + self.decimation_factor, + ) def forward_run(self): # set up monitors @@ -361,6 +380,7 @@ def adjoint_run(self): self.adjoint_design_region_monitors = [] self.adjoint_source_monitors = [] + self.adjoint_geometry_monitors = [] for ar in range(len(self.objective_functions)): # Reset the fields self.sim.restart_fields() @@ -381,6 +401,15 @@ def adjoint_run(self): # register a monitor over each differentiable source's support; the # adjoint field there is the gradient with respect to its currents + self.adjoint_geometry_monitors.append( + geometry_gradient.install_geometry_monitors( + self.sim, + self.differentiable_objects, + self.frequencies, + self.decimation_factor, + ) + ) + self.adjoint_source_monitors.append( source_gradient.install_source_gradient_monitors( self.sim, @@ -498,8 +527,31 @@ def _source_gradient_for(self, src, monitors, scale): ) return grads + def calculate_geometry_gradient(self): + """Shape derivatives for each flagged geometric object.""" + if not self.differentiable_objects: + return {} + out = {} + for ar in range(len(self.objective_functions)): + for oi, (object_index, obj) in enumerate(self.differentiable_objects): + grads = geometry_gradient.gradient( + self.sim, + obj, + object_index, + self.forward_geometry_monitors[oi], + self.adjoint_geometry_monitors[ar][oi], + self.frequencies, + ) + key = geometry_gradient.object_key(obj, object_index) + if len(self.objective_functions) == 1: + out[key] = grads + else: + out.setdefault(key, []).append(grads) + return out + def calculate_gradient(self): self.source_gradient = self.calculate_source_gradient() + self.geometry_gradient = self.calculate_geometry_gradient() # Iterate through all design regions and calculate gradient self.gradient = [ diff --git a/python/geom.py b/python/geom.py index b1207c75b..b15da1ce2 100755 --- a/python/geom.py +++ b/python/geom.py @@ -1033,6 +1033,35 @@ def __init__( self.pumping_rate = pumping_rate +# Parameters a geometric object can be differentiated with respect to. +# Unlike a source, where `center` and `size` move which grid points are driven +# and so fall outside the formulation, epsilon is a function of position and +# moving an object is exactly what a derivative with respect to position means. +_DIFFERENTIABLE_GEOMETRY = ("center", "size") + + +def _validate_differentiable_geometry(differentiable): + """Check and normalize a geometric object's `differentiable` argument.""" + if differentiable is None: + return () + if isinstance(differentiable, str): + raise ValueError( + "`differentiable` takes a list of parameter names, not a bare " + f"string; use ['{differentiable}'] instead." + ) + names = list(differentiable) + for name in names: + if name not in _DIFFERENTIABLE_GEOMETRY: + raise ValueError( + f"'{name}' is not a differentiable parameter of a geometric " + f"object. Valid choices are: " + f"{', '.join(_DIFFERENTIABLE_GEOMETRY)}." + ) + if len(set(names)) != len(names): + raise ValueError(f"`differentiable` contains duplicate names: {names}") + return tuple(names) + + class GeometricObject: """ This class, and its descendants, are used to specify the solid geometric objects that @@ -1088,7 +1117,13 @@ class GeometricObject: """ def __init__( - self, material=Medium(), center=Vector3(), epsilon_func=None, label=None + self, + material=Medium(), + center=Vector3(), + epsilon_func=None, + label=None, + differentiable=None, + name=None, ): """ Construct a `GeometricObject`. @@ -1119,6 +1154,8 @@ def __init__( self.label = label self.material = material + self.name = name + self.differentiable = _validate_differentiable_geometry(differentiable) self.center = Vector3(*center) def __contains__(self, point): diff --git a/python/meep.i b/python/meep.i index 51fa5aa3d..8b331eb97 100644 --- a/python/meep.i +++ b/python/meep.i @@ -652,6 +652,7 @@ void _get_eigenmode(meep::fields *f, double frequency, meep::direction d, const %feature("nothreadallow") py_do_harminv; %feature("nothreadallow") _get_array_slice_dimensions; %feature("nothreadallow") _get_gradient; +%feature("nothreadallow") _get_geometry_gradient; %feature("nothreadallow") _get_dft_array; %numpy_typemaps(std::complex, NPY_CDOUBLE, int); @@ -866,6 +867,47 @@ meep::volume_list *make_volume_list(const meep::volume &v, int c, // typemaps needed for material grid //-------------------------------------------------- +%inline %{ +void _get_geometry_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, int object_index, PyObject *params, + double fd_step) { + + PyArrayObject *pao_grad = (PyArrayObject *)grad; + if (!PyArray_Check(pao_grad)) meep::abort("grad parameter must be numpy array."); + if (!PyArray_ISCARRAY(pao_grad)) meep::abort("Numpy grad array must be C-style contiguous."); + if (PyArray_NDIM(pao_grad) != 2) meep::abort("Numpy grad array must have 2 dimensions."); + double *grad_c = (double *)PyArray_DATA(pao_grad); + npy_intp nparams = PyArray_DIMS(pao_grad)[1]; + + PyArrayObject *pao_params = (PyArrayObject *)params; + if (!PyArray_Check(pao_params)) meep::abort("params must be a numpy array."); + if (!PyArray_ISCARRAY(pao_params)) meep::abort("Numpy params array must be C-style contiguous."); + if (PyArray_DIMS(pao_params)[0] != nparams) + meep::abort("params has %td entries but grad is allocated for %td.", + PyArray_DIMS(pao_params)[0], nparams); + int *params_c = (int *)PyArray_DATA(pao_params); + + std::vector adjoint_fields = {fields_a_0,fields_a_1,fields_a_2}; + std::vector forward_fields = {fields_f_0,fields_f_1,fields_f_2}; + + PyArrayObject *pao_freqs = (PyArrayObject *)frequencies; + if (!PyArray_Check(pao_freqs)) meep::abort("frequencies parameter must be numpy array."); + if (!PyArray_ISCARRAY(pao_freqs)) meep::abort("Numpy frequencies array must be C-style contiguous."); + double *frequencies_c = (double *)PyArray_DATA(pao_freqs); + npy_intp nf = PyArray_DIMS(pao_freqs)[0]; + 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); + + meep_geom::geometry_addgradient(grad_c, nparams, nf, adjoint_fields, forward_fields, + frequencies_c, scalegrad, *grid_volume, geps, + object_index, params_c, fd_step); +} +%} + %inline %{ 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, diff --git a/python/tests/test_geometry_gradient.py b/python/tests/test_geometry_gradient.py new file mode 100644 index 000000000..3a3d24595 --- /dev/null +++ b/python/tests/test_geometry_gradient.py @@ -0,0 +1,264 @@ +"""Tests for adjoint gradients with respect to a geometric object's shape. + +Two conventions run through all of these, and both cost real debugging time to +establish. + +Object faces are placed at *quarter*-pixel offsets. Yee components sit half a +pixel apart, so a face at a pixel centre for one component lies exactly on a +pixel edge for another. On an edge no pixel straddles the face: the filling +fraction is 0 or 1 on both sides, meep applies no smoothing there, and the +derivative is one-sided. A quarter-pixel offset avoids both integer and +half-integer multiples of the pixel, so every component straddles every face +and the derivative is two-sided everywhere. + +Run lengths are pinned rather than left to `stop_when_dft_decayed`, for the +same reason as the source-gradient tests: an adaptive stop makes a perturbed +run end at a different time, and the difference scales with the perturbation. +""" + +import unittest +import warnings + +import numpy as np +from autograd import numpy as npa + +import meep as mp +import meep.adjoint as mpa + +RES = 20 +FCEN = 1.0 +RUN = 200.0 +CELL = mp.Vector3(8, 6) +SRC_C = mp.Vector3(-2.0, 0) +MON_C = mp.Vector3(2.0, 0) +RHO = 0.5 * np.ones(64) +FD_STEP = 1e-3 + + +def quarter_pixel(resolution: float) -> float: + """The offset that keeps a face clear of every component's pixel edges.""" + return 0.25 / resolution + + +def _design_region(sim): + """An inert design region; OptimizationProblem requires one.""" + grid = mp.MaterialGrid( + mp.Vector3(8, 8), mp.air, mp.Medium(index=1.5), grid_type="U_MEAN" + ) + region = mpa.DesignRegion( + grid, volume=mp.Volume(center=mp.Vector3(0, -2.0), size=mp.Vector3(0.4, 0.4)) + ) + sim.geometry = list(sim.geometry) + [ + mp.Block(center=region.center, size=region.size, material=grid) + ] + return region + + +def _problem(block, res=RES, component=mp.Hz, eps_averaging=True): + sim = mp.Simulation( + cell_size=CELL, + resolution=res, + boundary_layers=[mp.PML(1.0)], + sources=[ + mp.Source( + mp.GaussianSource(FCEN, fwidth=0.2), component=component, center=SRC_C + ) + ], + geometry=[block], + eps_averaging=eps_averaging, + force_complex_fields=True, + ) + region = _design_region(sim) + monitor = mpa.FourierFields( + sim, mp.Volume(center=MON_C, size=mp.Vector3(0, 0)), component + ) + return mpa.OptimizationProblem( + simulation=sim, + objective_functions=[lambda f: npa.sum(npa.abs(f) ** 2)], + objective_arguments=[monitor], + design_regions=[region], + frequencies=[FCEN], + minimum_run_time=RUN, + maximum_run_time=RUN, + ) + + +def _block(center, size, index=2.5, differentiable=None, res=RES): + """A block whose faces avoid every component's pixel edges.""" + offset = quarter_pixel(res) + return mp.Block( + center=mp.Vector3(center.x, center.y + offset), + size=mp.Vector3(size.x + 2 * offset, size.y), + material=mp.Medium(index=index), + differentiable=differentiable, + name="scatterer", + ) + + +class TestDifferentiableFlag(unittest.TestCase): + """Validation of `differentiable=` on a geometric object.""" + + def test_accepts_center_and_size(self): + block = mp.Block( + center=mp.Vector3(), + size=mp.Vector3(1, 1), + differentiable=["center", "size"], + name="b", + ) + self.assertEqual(block.differentiable, ("center", "size")) + self.assertEqual(block.name, "b") + + def test_default_is_not_differentiable(self): + self.assertEqual(mp.Block(size=mp.Vector3(1, 1)).differentiable, ()) + + def test_rejects_unknown_parameter(self): + with self.assertRaisesRegex(ValueError, "not a differentiable parameter"): + mp.Block(size=mp.Vector3(1, 1), differentiable=["waist"]) + + def test_rejects_bare_string(self): + with self.assertRaisesRegex(ValueError, "list of parameter names"): + mp.Block(size=mp.Vector3(1, 1), differentiable="center") + + def test_center_and_size_are_allowed_here_unlike_for_sources(self): + # A source's cotangent is gathered over a fixed set of grid points, so + # moving it changes which points it occupies rather than the amplitudes + # applied to them, and 'center' is rejected there. Epsilon is a + # function of position, so for geometry the same name is exactly what a + # derivative with respect to position means. + with self.assertRaises(ValueError): + mp.Source( + mp.GaussianSource(FCEN, fwidth=0.2), + component=mp.Ez, + center=SRC_C, + differentiable=["center"], + ) + mp.Block(size=mp.Vector3(1, 1), differentiable=["center"]) + + +class TestGuards(unittest.TestCase): + def test_refuses_without_subpixel_smoothing(self): + # Without smoothing the permittivity is a step function of position: + # nothing changes until a boundary crosses a pixel edge, then it + # changes by the full contrast. A finite difference of that is not a + # derivative, and returning a number would hide it. + block = _block(mp.Vector3(), mp.Vector3(1.0, 0.8), differentiable=["center"]) + opt = _problem(block, eps_averaging=False) + with self.assertRaisesRegex(ValueError, "subpixel smoothing"): + opt([RHO]) + + def test_warns_when_a_face_lies_on_a_pixel_edge(self): + # Pixel centres are at integer multiples of dx, so edges are at + # half-integers. A face there is straddled by nothing. + half = 0.5 / RES + block = mp.Block( + center=mp.Vector3(0, half), + size=mp.Vector3(1.0, 0.8), + material=mp.Medium(index=2.5), + differentiable=["center"], + name="s", + ) + opt = _problem(block) + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + opt([RHO]) + self.assertTrue( + any("pixel edge" in str(w.message) for w in caught), + "expected a warning about a face on a pixel edge", + ) + + +class TestShapeDerivative(unittest.TestCase): + SIZE = mp.Vector3(1.0, 0.8) + + def _value(self, center, size, res=RES, index=2.5): + opt = _problem(_block(center, size, index, res=res), res=res) + return float(np.asarray(opt([RHO], need_gradient=False)[0]).item()) + + def _adjoint(self, names, res=RES, index=2.5): + opt = _problem( + _block(mp.Vector3(), self.SIZE, index, differentiable=names, res=res), + res=res, + ) + _, grad = opt([RHO]) + return grad["scatterer"] + + def test_center_gradient(self): + for label, res, index in ( + ("baseline", RES, 2.5), + ("higher resolution", 40, 2.5), + ("lower contrast", RES, 1.5), + ): + with self.subTest(label): + gradient = self._adjoint(["center"], res, index)["center"] + adjoint = float(np.real(np.atleast_1d(gradient)[1])) + d = mp.Vector3(y=FD_STEP) + reference = ( + self._value(d, self.SIZE, res, index) + - self._value(mp.Vector3() - d, self.SIZE, res, index) + ) / (2 * FD_STEP) + self.assertLess( + abs(adjoint - reference) / abs(reference), + 2e-2, + f"{label}: adjoint {adjoint} vs finite difference {reference}", + ) + + def test_size_gradient(self): + adjoint = float(np.real(np.atleast_1d(self._adjoint(["size"])["size"])[1])) + d = mp.Vector3(y=FD_STEP) + reference = ( + self._value(mp.Vector3(), self.SIZE + d) + - self._value(mp.Vector3(), self.SIZE - d) + ) / (2 * FD_STEP) + self.assertLess(abs(adjoint - reference) / abs(reference), 2e-2) + + def test_translation_through_a_uniform_medium_is_zero(self): + # A block of the background material is not there at all, so moving it + # cannot change anything. Exact, and independent of finite-difference + # error, which makes it sharper than the comparisons above. + gradient = self._adjoint(["center"], index=1.0)["center"] + self.assertLess(float(np.max(np.abs(np.atleast_1d(gradient)))), 1e-12) + + def test_returns_one_entry_per_named_parameter(self): + gradient = self._adjoint(["center", "size"]) + self.assertEqual(set(gradient), {"center", "size"}) + for name in ("center", "size"): + self.assertEqual(np.shape(np.atleast_1d(gradient[name])), (3,)) + + def test_object_with_an_infinite_extent(self): + # `mp.inf` is the idiomatic way to write a layer that spans the cell, + # and it is 1e20 rather than a flag. A monitor built around one without + # clamping fails inside meep with "impossible(?) looping boundaries", + # which is what a stratified stack -- a metal reflector under a grating, + # say -- runs into immediately. + offset = quarter_pixel(RES) + + def slab(dy, differentiable=None): + return mp.Block( + center=mp.Vector3(0, dy + offset), + size=mp.Vector3(mp.inf, 0.4), + material=mp.Medium(index=2.5), + differentiable=differentiable, + name="slab", + ) + + opt = _problem(slab(0.0, differentiable=["center"])) + _, grad = opt([RHO]) + adjoint = float(np.real(np.atleast_1d(grad["slab"]["center"])[1])) + + def value(dy): + return float( + np.asarray(_problem(slab(dy))([RHO], need_gradient=False)[0]).item() + ) + + reference = (value(FD_STEP) - value(-FD_STEP)) / (2 * FD_STEP) + self.assertGreater(abs(reference), 1e-6, "objective must respond to the move") + self.assertLess(abs(adjoint - reference) / abs(reference), 2e-2) + + def test_no_flagged_object_leaves_the_return_shape_alone(self): + opt = _problem(_block(mp.Vector3(), self.SIZE)) + _, gradient = opt([RHO]) + self.assertNotIsInstance(gradient, dict) + + +if __name__ == "__main__": + unittest.main() diff --git a/src/meepgeom.cpp b/src/meepgeom.cpp index dce99195e..96c595350 100644 --- a/src/meepgeom.cpp +++ b/src/meepgeom.cpp @@ -1066,7 +1066,7 @@ void geom_epsilon::eff_chi1inv_row(meep::component c, double chi1inv_row[3], con void geom_epsilon::eff_chi1inv_matrix(meep::component c, symm_matrix *chi1inv_matrix, const meep::volume &v, double tol, int maxeval, - bool &fallback) { + bool &fallback, double fill_override) { const geometric_object *o; material_type mat, mat_behind; symm_matrix meps; @@ -1107,7 +1107,8 @@ void geom_epsilon::eff_chi1inv_matrix(meep::component c, symm_matrix *chi1inv_ma pixel.low = vector3_minus(pixel.low, shiftby); pixel.high = vector3_minus(pixel.high, shiftby); - double fill = box_overlap_with_object(pixel, *o, tol, maxeval); + double fill = + fill_override >= 0.0 ? fill_override : box_overlap_with_object(pixel, *o, tol, maxeval); material_epsmu(meep::type(c), mat, &meps, chi1inv_matrix); symm_matrix eps2, epsinv2; @@ -1198,6 +1199,29 @@ void geom_epsilon::eff_chi1inv_matrix(meep::component c, symm_matrix *chi1inv_ma sym_matrix_invert(chi1inv_matrix, &meps); } +/* The filling fraction this pixel would smooth with. Returns false when the + pixel does not straddle an interface between two distinct materials, in + which case the shape derivative vanishes there and the caller can skip it -- + which is what restricts the gradient loop to the boundary shell. */ +bool geom_epsilon::interface_fill(const meep::volume &v, double tol, int maxeval, double &fill, + const geometric_object **which, vector3 *shift) { + const geometric_object *o; + material_type mat, mat_behind; + vector3 p, shiftby; + fill = -1.0; + + if (!get_front_object(v, geometry_tree, p, &o, shiftby, mat, mat_behind)) return false; + if (material_type_equal(mat, mat_behind)) return false; + + geom_box pixel = gv2box(v); + pixel.low = vector3_minus(pixel.low, shiftby); + pixel.high = vector3_minus(pixel.high, shiftby); + fill = box_overlap_with_object(pixel, *o, tol, maxeval); + if (which) *which = o; + if (shift) *shift = shiftby; + return true; +} + static int eps_ever_negative = 0; static meep::field_type func_ft = meep::E_stuff; @@ -2940,6 +2964,315 @@ static std::complex forward_dft_value(const meep::dft_chunk *ch, if (i >= 0 && (size_t)i < ch->N) return ch->dft[nf * i + f_i]; } return 0; + +/* ------------------------------------------------------------------ */ +/* Gradients with respect to a geometric object's centre and size. */ +/* ------------------------------------------------------------------ */ + +/* The smoothed permittivity depends on the geometry only through two things: + the filling fraction of each pixel and the interface normal. So + + d(chi1inv)/d(parameter) = d(chi1inv)/d(fill) * d(fill)/d(parameter) + + with the second factor analytic for an axis-aligned block, since the block + is an intersection of three slabs and the overlap volume factorizes: + + |pixel ∩ block| = prod_i overlap_i(c_i, s_i) + overlap_i = clamp(min(c_i + s_i/2, hi_i) - max(c_i - s_i/2, lo_i), 0, dx) + + d(overlap_i)/d(c_i) is -1, 0 or +1 depending on which face lies inside the + pixel, and d(overlap_i)/d(s_i) is 0 or +/-1/2. Everything else is a product + of the other axes' overlaps. + + The first factor comes from varying `fill` through eff_chi1inv_matrix rather + than re-deriving Kottke's algebra here, which keeps this consistent with + whatever meep actually does. That variation is pure local algebra: `delta` + is linear in fill by construction, and only the final delta -> chi1inv map + is nonlinear. + + Perturbing the geometry and re-differencing instead -- which is what this + routine used to do -- means differencing box_overlap_with_object, an + adaptive quadrature, so its tolerance is amplified by 1/step. That produced + a gradient with no stable step regime at all. + + Only pixels straddling the boundary contribute, since d(fill)/d(parameter) + vanishes wherever a pixel is wholly inside or wholly outside. That is the + discrete form of a shape derivative being a surface integral, and it is why + the loop below skips interior pixels. */ + +enum geom_param { + GEOM_CENTER_X = 0, + GEOM_CENTER_Y, + GEOM_CENTER_Z, + GEOM_SIZE_X, + GEOM_SIZE_Y, + GEOM_SIZE_Z +}; + +/* Overlap of [lo, hi] with the block's extent along one axis, and its + derivatives with respect to that axis's centre and half-size. */ +static void axis_overlap(double lo, double hi, double c, double s, double &overlap, + double &d_dcenter, double &d_dsize) { + const double blo = c - 0.5 * s, bhi = c + 0.5 * s; + const double left = std::max(lo, blo), right = std::min(hi, bhi); + overlap = right - left; + if (overlap <= 0) { + overlap = 0; + d_dcenter = d_dsize = 0; + return; + } + /* Moving the centre moves both faces together; growing the size moves them + apart by half each. A face contributes only in the pixel that contains it. + + The comparisons are half-open on purpose. With strict inequalities on both + sides, a face lying exactly on a pixel boundary belongs to neither + neighbour and its contribution vanishes -- silently, and for every pixel + along that face, which zeroes the whole derivative for that axis. Round + geometry on a round grid hits this constantly: a block of width 1.0 at + resolution 20 has faces exactly 10 pixels from its centre. + + Half-open assigns such a face to exactly one of the two pixels, so nothing + is dropped and nothing is double counted. The derivative is then one-sided + there, which is the truth: the smoothed permittivity has a kink at a pixel + boundary and no two-sided derivative exists. */ + const double lower_inside = (blo >= lo && blo < hi) ? 1.0 : 0.0; + const double upper_inside = (bhi > lo && bhi <= hi) ? 1.0 : 0.0; + d_dcenter = upper_inside - lower_inside; + d_dsize = 0.5 * (upper_inside + lower_inside); +} + +void geometry_addgradient(double *v, size_t nparams, size_t nf, + std::vector fields_a, + std::vector fields_f, double *frequencies, + double scalegrad, meep::grid_volume &gv, geom_epsilon *geps, + int object_index, int *params, double du) { + (void)frequencies; + (void)du; + if (object_index < 0 || object_index >= geps->geometry.num_items) + meep::abort("geometry_addgradient: object index %d out of range (%d objects)", object_index, + geps->geometry.num_items); + geometric_object *obj = &geps->geometry.items[object_index]; + if (obj->which_subclass != geometric_object::BLOCK) + meep::abort("geometry_addgradient: only blocks are supported"); + + const vector3 bc = obj->center; + const vector3 bs = obj->subclass.block_data->size; + const double centers[3] = {bc.x, bc.y, bc.z}; + const double sizes[3] = {bs.x, bs.y, bs.z}; + + std::vector adjoint_chunks[3], forward_chunks[3]; + for (int i = 0; i < 3; i++) { + for (meep::dft_chunk *c = fields_a[i]->chunks; c; c = c->next_in_dft) + adjoint_chunks[i].push_back(c); + for (meep::dft_chunk *c = fields_f[i]->chunks; c; c = c->next_in_dft) + forward_chunks[i].push_back(c); + } + + std::vector local(nf * nparams, 0.0); + /* `fill` is dimensionless and O(1) and `delta` is linear in it, so this step + is well conditioned -- unlike a step in a length, which has to be compared + against the pixel. */ + const double dfill = 1e-6; + + for (size_t f_i = 0; f_i < nf; f_i++) { + for (int ci_adjoint = 0; ci_adjoint < 3; ci_adjoint++) { + int num_chunks = adjoint_chunks[ci_adjoint].size(); + if (num_chunks == 0) continue; + + for (int cur = 0; cur < num_chunks; cur++) { + meep::dft_chunk *adj_chunk = adjoint_chunks[ci_adjoint][cur]; + meep::component adjoint_c = adj_chunk->c; + meep::grid_volume gv_adj = gv.subvolume(adj_chunk->is, adj_chunk->ie, adjoint_c); + + for (int ci_forward = 0; ci_forward < 3; ci_forward++) { + /* Pair the forward chunk with this adjoint chunk *spatially*. + Indexing both lists by `cur` assumes they are the same length and + in the same order, and they need not be -- which is why + matching_dft_chunk exists a few functions up, for exactly this + purpose in material_grids_addgradient. A foreign chunk here does + not merely mix the gradient up: it then gets indexed with this + chunk's offsets and read past the end of its array. */ + if (forward_chunks[ci_forward].empty()) continue; + meep::dft_chunk *fwd_chunk = + matching_dft_chunk(forward_chunks[ci_forward], adj_chunk); + if (!fwd_chunk) continue; + meep::component forward_c = fwd_chunk->c; + meep::grid_volume gv_fwd = gv.subvolume(fwd_chunk->is, fwd_chunk->ie, forward_c); + + int dir_idx; + switch (meep::component_direction(forward_c)) { + case meep::X: + case meep::R: dir_idx = 0; break; + case meep::Y: + case meep::P: dir_idx = 1; break; + case meep::Z: dir_idx = 2; break; + default: continue; + } + + meep::ivec loop_is = adj_chunk->persist ? adj_chunk->is_old : adj_chunk->is; + meep::ivec loop_ie = adj_chunk->persist ? adj_chunk->ie_old : adj_chunk->ie; + + LOOP_OVER_IVECS(gv_adj, loop_is, loop_ie, idx_adj) { + IVEC_LOOP_ILOC(gv_adj, ip); + IVEC_LOOP_LOC(gv_adj, p); + std::complex adj = adj_chunk->dft[nf * idx_adj + f_i]; + if (adj == 0.0) continue; + + /* The smoothed tensor is diagonal in Cartesian axes only where the + interface normal lies along an axis, which for a block is true on + a face and false on an edge or corner. Those pixels are O(N) of + the O(N^2) on the boundary, so their share falls with resolution + -- which is the signature of the residual error here. Contracting + them needs the forward field of another component, restricted to + the two epsilon nodes between the pair. */ + std::complex fwd; + meep::vec eps_at = p; + double node_weight = 1.0; + int num_nodes = 1; + meep::vec node_pos[2]; + std::complex node_fwd[2]; + if (forward_c == adjoint_c) { + node_pos[0] = p; + /* `idx_adj` is an offset into gv_adj, so it cannot be used to + index the forward chunk. Look the value up by location, which + bounds-checks it as the two lookups below already did. */ + node_fwd[0] = + forward_dft_value(fwd_chunk, fwd_chunk->dft, gv_fwd, gv, ip, nf, f_i); + } + else { + num_nodes = 2; + node_weight = 0.5; + meep::ivec fwd_p = ip + gv.iyee_shift(forward_c) - gv.iyee_shift(adjoint_c); + meep::ivec unit_a = unit_ivec(gv.dim, component_direction(adjoint_c)); + meep::ivec unit_f = unit_ivec(gv.dim, component_direction(forward_c)); + meep::ivec pl[2] = {fwd_p, fwd_p + unit_a * 2}; + meep::ivec pr[2] = {fwd_p - unit_f * 2, fwd_p + unit_a * 2 - unit_f * 2}; + for (int nd = 0; nd < 2; nd++) { + ptrdiff_t i1 = gv_fwd.index(forward_c, pl[nd]); + ptrdiff_t i2 = gv_fwd.index(forward_c, pr[nd]); + std::complex f1 = + ((i1 >= fwd_chunk->N) || (i1 < 0)) ? 0 : fwd_chunk->dft[nf * i1 + f_i]; + std::complex f2 = + ((i2 >= fwd_chunk->N) || (i2 < 0)) ? 0 : fwd_chunk->dft[nf * i2 + f_i]; + node_fwd[nd] = std::complex(0.5, 0) * (f1 + f2); + node_pos[nd] = gv[(pl[nd] + pr[nd]) / 2]; + } + } + + for (int node = 0; node < num_nodes; node++) { + fwd = node_fwd[node]; + if (fwd == 0.0) continue; + const meep::vec p_node = node_pos[node]; + + meep::volume voxel(p_node); + LOOP_OVER_DIRECTIONS(gv.dim, d) { + voxel.set_direction_min(d, p_node.in_direction(d) - 0.5 * gv.inva); + voxel.set_direction_max(d, p_node.in_direction(d) + 0.5 * gv.inva); + } + + /* Skip pixels that do not straddle the boundary: d(fill)/d(param) + is zero there, so they contribute nothing. */ + /* `fill` is the fraction of the pixel inside whichever object + get_front_object selected. If that is not the object being + differentiated, d(fill)/d(our parameters) is not what this pixel + responds to -- and taking it anyway gets the sign wrong wherever + the front object is the background instead. */ + double fill; + const geometric_object *front = NULL; + vector3 shiftby = {0, 0, 0}; + if (!geps->interface_fill(voxel, geps->tol, geps->maxeval, fill, &front, &shiftby)) + continue; + if (front != obj) continue; + if (fill <= 0.0 || fill >= 1.0) continue; + + /* d(fill)/d(parameter), analytic and separable. */ + double overlap[3], d_dc[3], d_ds[3], extent[3]; + bool degenerate = false; + for (int ax = 0; ax < 3; ax++) { + meep::direction dd = (ax == 0) ? meep::X : (ax == 1) ? meep::Y : meep::Z; + const bool resolved = (gv.dim == meep::D3) || (gv.dim == meep::D2 && ax < 2) || + (gv.dim == meep::D1 && ax == 2); + if (!resolved) { + /* A dimension the simulation does not resolve. The block is + effectively infinite along it, so it contributes a factor of + one to the overlap and nothing to the derivative -- note the + block's own size along such an axis is typically zero, so + using it here would make every pixel degenerate. */ + extent[ax] = 1.0; + overlap[ax] = 1.0; + d_dc[ax] = d_ds[ax] = 0.0; + continue; + } + const double sh = (ax == 0) ? shiftby.x : (ax == 1) ? shiftby.y : shiftby.z; + const double lo = voxel.in_direction_min(dd) - sh; + const double hi = voxel.in_direction_max(dd) - sh; + extent[ax] = hi - lo; + axis_overlap(lo, hi, centers[ax], sizes[ax], overlap[ax], d_dc[ax], d_ds[ax]); + if (overlap[ax] <= 0) degenerate = true; + } + if (degenerate) continue; + + double pixel_volume = 1.0; + for (int ax = 0; ax < 3; ax++) + pixel_volume *= extent[ax]; + + /* d(chi1inv)/d(fill), from meep's own tensor assembly. */ + bool fb_lo = false, fb_hi = false; + symm_matrix m_lo, m_hi; + geps->eff_chi1inv_matrix(adjoint_c, &m_lo, voxel, geps->tol, geps->maxeval, fb_lo, + std::max(0.0, fill - dfill)); + geps->eff_chi1inv_matrix(adjoint_c, &m_hi, voxel, geps->tol, geps->maxeval, fb_hi, + std::min(1.0, fill + dfill)); + if (fb_lo || fb_hi) continue; + + const double lo_row[3] = {m_lo.m00, m_lo.m01, m_lo.m02}; + const double hi_row[3] = {m_hi.m00, m_hi.m01, m_hi.m02}; + const double lo_row1[3] = {m_lo.m01, m_lo.m11, m_lo.m12}; + const double hi_row1[3] = {m_hi.m01, m_hi.m11, m_hi.m12}; + const double lo_row2[3] = {m_lo.m02, m_lo.m12, m_lo.m22}; + const double hi_row2[3] = {m_hi.m02, m_hi.m12, m_hi.m22}; + int row_of = 0; + switch (meep::component_direction(adjoint_c)) { + case meep::X: + case meep::R: row_of = 0; break; + case meep::Y: + case meep::P: row_of = 1; break; + case meep::Z: row_of = 2; break; + default: continue; + } + const double *lo_r = (row_of == 0) ? lo_row : (row_of == 1) ? lo_row1 : lo_row2; + const double *hi_r = (row_of == 0) ? hi_row : (row_of == 1) ? hi_row1 : hi_row2; + const double actual_dfill = std::min(1.0, fill + dfill) - std::max(0.0, fill - dfill); + if (actual_dfill <= 0) continue; + const double dchi_dfill = (hi_r[dir_idx] - lo_r[dir_idx]) / actual_dfill; + + const std::complex pair = + std::complex(double(adj.real()), double(adj.imag())) * + std::complex(double(fwd.real()), double(fwd.imag())); + const double cyl_scale = (gv.dim == meep::Dcyl) ? 2 * p_node.r() : 1; + + for (size_t ip = 0; ip < nparams; ip++) { + const int which = params[ip]; + const int ax = which % 3; + double dfill_dp = (which < 3) ? d_dc[ax] : d_ds[ax]; + if (dfill_dp == 0.0) continue; + for (int other = 0; other < 3; other++) + if (other != ax) dfill_dp *= overlap[other]; + dfill_dp /= pixel_volume; + + /* the leading minus matches get_material_gradient's convention, + which returns -(d row/d parameter) */ + local[nparams * f_i + ip] -= + node_weight * scalegrad * cyl_scale * dchi_dfill * dfill_dp * std::real(pair); + } + } // node + } + } + } + } + } + + meep::sum_to_all(local.data(), v, int(nf * nparams)); } void material_grids_addgradient(double *v, size_t ng, size_t nf, diff --git a/src/meepgeom.hpp b/src/meepgeom.hpp index 004954209..dee643ec8 100644 --- a/src/meepgeom.hpp +++ b/src/meepgeom.hpp @@ -212,8 +212,23 @@ class geom_epsilon : public meep::material_function { virtual bool is_thread_safe() const { return !has_user_materials; } bool has_user_materials; + /* `fill_override >= 0` substitutes that filling fraction for the one this + routine would compute from the geometry. The shape derivative uses it: the + smoothed tensor depends on the geometry only through `fill` and the + interface normal, so d(chi1inv)/d(parameter) factors into + d(chi1inv)/d(fill) -- obtained by varying `fill` here, which is pure local + algebra with no quadrature -- times d(fill)/d(parameter), which is + analytic for an axis-aligned block. Differencing the geometry directly + instead means differencing box_overlap_with_object, an adaptive + quadrature, which amplifies its tolerance by 1/step. */ void eff_chi1inv_matrix(meep::component c, symm_matrix *chi1inv_matrix, const meep::volume &v, - double tol, int maxeval, bool &fallback); + double tol, int maxeval, bool &fallback, double fill_override = -1.0); + + /* The filling fraction and interface normal this pixel would smooth with, + or fill < 0 if the pixel does not straddle an interface (so the shape + derivative is zero there and the pixel can be skipped). */ + bool interface_fill(const meep::volume &v, double tol, int maxeval, double &fill, + const geometric_object **which = NULL, vector3 *shift = NULL); void fallback_chi1inv_row(meep::component c, double chi1inv_row[3], const meep::volume &v, double tol, int maxeval); @@ -294,6 +309,12 @@ meep::vec material_grid_grad(vector3 p, material_data *md, const geometric_objec double matgrid_val(vector3 p, geom_box_tree tp, int oi, material_data *md); double material_grid_val(vector3 p, material_data *md); geom_box_tree calculate_tree(const meep::volume &v, geometric_object_list g); +void geometry_addgradient(double *v, size_t nparams, size_t nf, + std::vector fields_a, + std::vector fields_f, double *frequencies, + double scalegrad, meep::grid_volume &gv, geom_epsilon *geps, + int object_index, int *params, double du); + void material_grids_addgradient(double *v, size_t ng, size_t nf, std::vector fields_a, std::vector fields_f, double *frequencies,