diff --git a/python/adjoint/__init__.py b/python/adjoint/__init__.py index 7a898e677..949d3d689 100644 --- a/python/adjoint/__init__.py +++ b/python/adjoint/__init__.py @@ -5,9 +5,12 @@ from .objective import * +from . import utils +from .utils import DesignRegion + from .basis import BilinearInterpolationBasis -from .optimization_problem import (OptimizationProblem, Grid, DesignRegion) +from .optimization_problem import (OptimizationProblem) from .filter_source import FilteredSource @@ -15,8 +18,6 @@ from .connectivity import * -from . import utils - try: from .wrapper import MeepJaxWrapper except ModuleNotFoundError as _: diff --git a/python/adjoint/objective.py b/python/adjoint/objective.py index 11d78e778..9789d94b8 100644 --- a/python/adjoint/objective.py +++ b/python/adjoint/objective.py @@ -4,9 +4,10 @@ import numpy as np import meep as mp from .filter_source import FilteredSource -from .optimization_problem import Grid from meep.simulation import py_v3_to_vec +from collections import namedtuple +Grid = namedtuple('Grid', ['x', 'y', 'z', 'w']) class ObjectiveQuantity(abc.ABC): """A differentiable objective quantity. diff --git a/python/adjoint/optimization_problem.py b/python/adjoint/optimization_problem.py index c33ecb094..b3d0577be 100644 --- a/python/adjoint/optimization_problem.py +++ b/python/adjoint/optimization_problem.py @@ -3,48 +3,7 @@ from autograd import grad, jacobian from collections import namedtuple -Grid = namedtuple('Grid', ['x', 'y', 'z', 'w']) -YeeDims = namedtuple('YeeDims', ['Ex', 'Ey', 'Ez']) - - -class DesignRegion(object): - def __init__( - self, - design_parameters, - volume=None, - size=None, - center=mp.Vector3(), - MaterialGrid=None, - ): - self.volume = volume if volume else mp.Volume(center=center, size=size) - self.size = self.volume.size - self.center = self.volume.center - self.design_parameters = design_parameters - self.num_design_params = design_parameters.num_params - self.MaterialGrid = MaterialGrid - - def update_design_parameters(self, design_parameters): - self.design_parameters.update_weights(design_parameters) - - def get_gradient(self, sim, fields_a, fields_f, frequencies): - for c in range(3): - fields_a[c] = fields_a[c].flatten(order='C') - fields_f[c] = fields_f[c].flatten(order='C') - fields_a = np.concatenate(fields_a) - fields_f = np.concatenate(fields_f) - num_freqs = np.array(frequencies).size - - grad = np.zeros((num_freqs, self.num_design_params)) # preallocate - - geom_list = sim.geometry - f = sim.fields - vol = sim._fit_volume_to_simulation(self.volume) - # compute the gradient - sim_is_cylindrical = (sim.dimensions == mp.CYLINDRICAL) or sim.is_cylindrical - mp._get_gradient(grad,fields_a,fields_f,vol,np.array(frequencies),geom_list,f, sim_is_cylindrical) - - return np.squeeze(grad).T - +from . import utils, DesignRegion class OptimizationProblem(object): """Top-level class in the MEEP adjoint module. @@ -76,9 +35,6 @@ def __init__( ): self.sim = simulation - self.components = [mp.Ex,mp.Ey,mp.Ez] - if self.sim.is_cylindrical or self.sim.dimensions == mp.CYLINDRICAL: - self.components = [mp.Er,mp.Ep,mp.Ez] if isinstance(objective_functions, list): self.objective_functions = objective_functions @@ -161,13 +117,13 @@ def __call__(self, rho_vector=None, need_value=True, need_gradient=True): print("Starting forward run...") self.forward_run() print("Starting adjoint run...") - self.a_E = [] + self.D_a = [] self.adjoint_run() print("Calculating gradient...") self.calculate_gradient() elif self.current_state == "FWD": print("Starting adjoint run...") - self.a_E = [] + self.D_a = [] self.adjoint_run() print("Calculating gradient...") self.calculate_gradient() @@ -210,28 +166,9 @@ def prepare_forward_run(self): self.forward_monitors.append(m.register_monitors(self.frequencies)) # register design region - self.design_region_monitors = [ - self.sim.add_dft_fields( - self.components, - self.frequencies, - where=dr.volume, - yee_grid=True, - decimation_factor=self.decimation_factor, - ) for dr in self.design_regions - ] - - # store design region voxel parameters - self.design_grids = [] - if self.sim.is_cylindrical or self.sim.dimensions == mp.CYLINDRICAL: - YeeDims = namedtuple('YeeDims', ['Er','Ep','Ez']) - else: - YeeDims = namedtuple('YeeDims', ['Ex','Ey','Ez']) - for drm in self.design_region_monitors: - s = [ - self.sim.get_array_slice_dimensions(c, vol=drm.where)[0] - for c in self.components - ] - self.design_grids += [YeeDims(*s)] + self.design_region_monitors = utils.install_design_region_monitors( + self.sim, self.design_regions, self.frequencies, self.decimation_factor + ) def forward_run(self): # set up monitors @@ -254,16 +191,8 @@ def forward_run(self): if len(self.f0) == 1: self.f0 = self.f0[0] - # Store forward fields for each set of design variables in array (x,y,z,field_components,frequencies) - self.d_E = [[ - np.zeros((self.nf, c[0], c[1], c[2]), dtype=np.complex128) - for c in dg - ] for dg in self.design_grids] - for nb, dgm in enumerate(self.design_region_monitors): - for ic, c in enumerate(self.components): - for f in range(self.nf): - self.d_E[nb][ic][f, :, :, :] = atleast_3d( - self.sim.get_dft_array(dgm, c, f)) + # Store forward fields for each set of design variables in array + self.D_f = utils.gather_design_region_fields(self.sim,self.design_region_monitors,self.frequencies) # store objective function evaluation in memory self.f_bank.append(self.f0) @@ -288,11 +217,14 @@ def adjoint_run(self): # set up adjoint sources and monitors self.prepare_adjoint_run() - if self.sim.is_cylindrical or self.sim.dimensions == mp.CYLINDRICAL: + # flip the m number + if utils._check_if_cylindrical(self.sim): self.sim.m = -self.sim.m + # flip the k point if self.sim.k_point: self.sim.k_point *= -1 + for ar in range(len(self.objective_functions)): # Reset the fields self.sim.reset_meep() @@ -301,15 +233,9 @@ def adjoint_run(self): self.sim.change_sources(self.adjoint_sources[ar]) # register design flux - self.design_region_monitors = [ - self.sim.add_dft_fields( - self.components, - self.frequencies, - where=dr.volume, - yee_grid=True, - decimation_factor=self.decimation_factor, - ) for dr in self.design_regions - ] + self.design_region_monitors = utils.install_design_region_monitors( + self.sim, self.design_regions, self.frequencies, self.decimation_factor + ) # Adjoint run self.sim.run(until_after_sources=mp.stop_when_dft_decayed( @@ -318,25 +244,17 @@ def adjoint_run(self): self.maximum_run_time )) - # Store adjoint fields for each design set of design variables in array (x,y,z,field_components,frequencies) - self.a_E.append([[ - np.zeros((self.nf, c[0], c[1], c[2]), dtype=np.complex128) - for c in dg - ] for dg in self.design_grids]) - for nb, dgm in enumerate(self.design_region_monitors): - for ic, c in enumerate(self.components): - for f in range(self.nf): - if (self.sim.is_cylindrical or self.sim.dimensions == mp.CYLINDRICAL): - # Addtional factor of 2 for cyldrical coordinate - self.a_E[ar][nb][ic][f, :, :, :] = 2 * atleast_3d(self.sim.get_dft_array(dgm, c, f)) - else: - self.a_E[ar][nb][ic][f, :, :, :] = atleast_3d(self.sim.get_dft_array(dgm, c, f)) - - if self.sim.is_cylindrical or self.sim.dimensions == mp.CYLINDRICAL: + # Store adjoint fields for each design set of design variables + self.D_a.append(utils.gather_design_region_fields(self.sim,self.design_region_monitors,self.frequencies)) + + # reset the m number + if utils._check_if_cylindrical(self.sim): self.sim.m = -self.sim.m - # update optimizer's state + # reset the k point if self.sim.k_point: self.sim.k_point *= -1 + + # update optimizer's state self.current_state = "ADJ" def calculate_gradient(self): @@ -344,8 +262,8 @@ def calculate_gradient(self): self.gradient = [[ dr.get_gradient( self.sim, - self.a_E[ar][dri], - self.d_E[dri], + self.D_a[ar][dri], + self.D_f[dri], self.frequencies, ) for dri, dr in enumerate(self.design_regions) ] for ar in range(len(self.objective_functions))] @@ -486,7 +404,7 @@ def calculate_fd_gradient( return fd_gradient, fd_gradient_idx - def update_design(self, rho_vector): + def update_design(self, rho_vector, beta=None): """Update the design permittivity function. rho_vector ....... a list of numpy arrays that maps to each design region @@ -497,6 +415,8 @@ def update_design(self, rho_vector): "Each vector of design variables must contain only one dimension." ) b.update_design_parameters(rho_vector[bi]) + if beta: + b.update_beta(beta) self.sim.reset_meep() self.current_state = "INIT" diff --git a/python/adjoint/utils.py b/python/adjoint/utils.py index decbe9e11..d5b5af0ab 100644 --- a/python/adjoint/utils.py +++ b/python/adjoint/utils.py @@ -3,14 +3,93 @@ import meep as mp import numpy as onp -from . import ObjectiveQuantity, DesignRegion +from . import ObjectiveQuantity # Meep field components used to compute adjoint sensitivities -_ADJOINT_FIELD_COMPONENTS = [mp.Ex, mp.Ey, mp.Ez] +_ADJOINT_FIELD_COMPONENTS = [mp.Dx, mp.Dy, mp.Dz] +_ADJOINT_FIELD_COMPONENTS_CYL = [mp.Dr, mp.Dp, mp.Dz] # The frequency axis in the array returned by `mp._get_gradient()` _GRADIENT_FREQ_AXIS = 1 +# The returned axis order from get_dft_array in cylindrical coordinates +_FREQ_AXIS = 0 +_RHO_AXIS = 2 +_PHI_AXIS = 3 +_Z_AXIS = 1 + +class DesignRegion(object): + def __init__( + self, + design_parameters, + volume=None, + size=None, + center=mp.Vector3() + ): + self.volume = volume if volume else mp.Volume(center=center, size=size) + self.size = self.volume.size + self.center = self.volume.center + self.design_parameters = design_parameters + self.num_design_params = design_parameters.num_params + + def update_design_parameters(self, design_parameters): + self.design_parameters.update_weights(design_parameters) + + def update_beta(self,beta): + self.design_parameters.beta=beta + + def get_gradient(self, sim, fields_a, fields_f, frequencies): + num_freqs = onp.array(frequencies).size + shapes = [] + '''We have the option to linear scale the gradients up front + using the scalegrad parameter (leftover from MPB API). Not + currently needed for any existing feature (but available for + future use)''' + scalegrad = 1 + for component_idx, component in enumerate(_compute_components(sim)): + '''We need to correct for the rare cases that get_dft_array + returns a singleton element for the forward and adjoint fields. + This only occurs when we are in 2D and only working with a particular + polarization (as the other fields are never stored). For example, the + 2D in-plane polarization consists of a single scalar Ez field + (technically, meep doesn't store anything for these cases, but get_dft_array + still returns 0). + + Our get_gradient algorithm, however, requires we pass an array of + zeros with the proper shape as the design_region.''' + spatial_shape = sim.get_array_slice_dimensions(component, vol=self.volume)[0] + if (fields_a[component_idx][0,...].size == 1): + fields_a[component_idx] = onp.zeros(onp.insert(spatial_shape,0,num_freqs)) + fields_f[component_idx] = onp.zeros(onp.insert(spatial_shape,0,num_freqs)) + if _check_if_cylindrical(sim): + '''For some reason, get_dft_array returns the field + components in a different order than the convention used + throughout meep. So, we reorder them here so we can use + the same field macros later in our get_gradient function. + ''' + fields_a[component_idx] = onp.transpose(fields_a[component_idx],(_FREQ_AXIS,_RHO_AXIS,_PHI_AXIS,_Z_AXIS)) + fields_f[component_idx] = onp.transpose(fields_f[component_idx],(_FREQ_AXIS,_RHO_AXIS,_PHI_AXIS,_Z_AXIS)) + shapes.append(fields_a[component_idx].shape) + fields_a[component_idx] = fields_a[component_idx].flatten(order='C') + fields_f[component_idx] = fields_f[component_idx].flatten(order='C') + shapes = onp.asarray(shapes).flatten(order='C') + fields_a = onp.concatenate(fields_a) + fields_f = onp.concatenate(fields_f) + + grad = onp.zeros((num_freqs, self.num_design_params)) # preallocate + geom_list = sim.geometry + f = sim.fields + vol = sim._fit_volume_to_simulation(self.volume) + + # compute the gradient + mp._get_gradient(grad,scalegrad,fields_a,fields_f,sim.gv,vol.swigobj,onp.array(frequencies),sim.geps,shapes) + return onp.squeeze(grad).T + +def _check_if_cylindrical(sim): + return sim.is_cylindrical or (sim.dimensions == mp.CYLINDRICAL) + +def _compute_components(sim): + return _ADJOINT_FIELD_COMPONENTS_CYL if _check_if_cylindrical(sim) else _ADJOINT_FIELD_COMPONENTS def _make_at_least_nd(x: onp.ndarray, dims: int = 3) -> onp.ndarray: """Makes an array have at least the specified number of dimensions.""" @@ -61,15 +140,16 @@ def install_design_region_monitors( simulation: mp.Simulation, design_regions: List[DesignRegion], frequencies: List[float], + decimation_factor: int = 0 ) -> List[mp.DftFields]: """Installs DFT field monitors at the design regions of the simulation.""" design_region_monitors = [ simulation.add_dft_fields( - _ADJOINT_FIELD_COMPONENTS, + _compute_components(simulation), frequencies, where=design_region.volume, yee_grid=True, - decimation_factor=0 + decimation_factor=decimation_factor ) for design_region in design_regions ] return design_region_monitors @@ -120,7 +200,7 @@ def gather_design_region_fields( fwd_fields = [] for monitor in design_region_monitors: fields_by_component = [] - for component in _ADJOINT_FIELD_COMPONENTS: + for component in _compute_components(simulation): fields_by_freq = [] for freq_idx, _ in enumerate(frequencies): fields = simulation.get_dft_array(monitor, component, freq_idx) diff --git a/python/adjoint/wrapper.py b/python/adjoint/wrapper.py index d36d1b89d..f918feeec 100644 --- a/python/adjoint/wrapper.py +++ b/python/adjoint/wrapper.py @@ -99,9 +99,7 @@ def __init__( monitors: List[EigenmodeCoefficient], design_regions: List[DesignRegion], frequencies: List[float], - measurement_interval: float = 50.0, - dft_field_components: Tuple[int, ...] = (mp.Ez, ), - dft_threshold: float = 1e-6, + dft_threshold: float = 1e-11, minimum_run_time: float = 0, maximum_run_time: float = onp.inf, until_after_sources: bool = True, @@ -111,14 +109,11 @@ def __init__( self.monitors = monitors self.design_regions = design_regions self.frequencies = frequencies - self.measurement_interval = measurement_interval - self.dft_field_components = dft_field_components self.dft_threshold = dft_threshold self.minimum_run_time = minimum_run_time self.maximum_run_time = maximum_run_time self.until_after_sources = until_after_sources - self._reset_convergence_measurement() self._simulate_fn = self._initialize_callable() def __call__(self, designs: List[jnp.ndarray]) -> jnp.ndarray: @@ -140,17 +135,6 @@ def __call__(self, designs: List[jnp.ndarray]) -> jnp.ndarray: """ return self._simulate_fn(designs) - def _reset_convergence_measurement(self, - monitors: List[mp.DftFields] = None - ) -> None: - """Resets the DFT convergence measurement.""" - if monitors is None: - monitors = [] - self._dft_convergence_monitors = monitors - self._last_measurement_meep_time = 0.0 - self._previous_fields = self._init_empty_dft_field_container() - self._dft_relative_change = [] - def _run_fwd_simulation(self, design_variables): """Runs forward simulation, returning monitor values and design region fields.""" utils.validate_and_update_design(self.design_regions, design_variables) @@ -165,9 +149,8 @@ def _run_fwd_simulation(self, design_variables): self.simulation.init_sim() sim_run_args = { 'until_after_sources' if self.until_after_sources else 'until': - self._simulation_run_callback + mp.stop_when_dft_decayed(self.dft_threshold,self.minimum_run_time,self.maximum_run_time) } - self._reset_convergence_measurement(design_region_monitors) self.simulation.run(**sim_run_args) monitor_values = utils.gather_monitor_values(self.monitors) @@ -197,9 +180,8 @@ def _run_adjoint_simulation(self, monitor_values_grad): self.simulation.init_sim() sim_run_args = { 'until_after_sources' if self.until_after_sources else 'until': - self._simulation_run_callback + mp.stop_when_dft_decayed(self.dft_threshold,self.minimum_run_time,self.maximum_run_time) } - self._reset_convergence_measurement(design_region_monitors) self.simulation.run(**sim_run_args) return utils.gather_design_region_fields( @@ -254,79 +236,3 @@ def _simulate_rev(res, monitor_values_grad): simulate.defvjp(_simulate_fwd, _simulate_rev) return simulate - - def _init_empty_dft_field_container(self) -> List[List[float]]: - """Initializes a nested list for storing DFT fields for convergence monitoring.""" - num_components = len(self.dft_field_components) - num_monitors = len(self._dft_convergence_monitors) - return [[0.0 for _ in range(num_components)] - for _ in range(num_monitors)] - - def _are_dfts_converged(self, sim: mp.Simulation) -> bool: - """Callback to determine whether the DFT fields are converged below the threshold.""" - if self.dft_threshold == 0 or not self._dft_convergence_monitors or not self.dft_field_components: - return False - relative_change = [] - current_fields = self._init_empty_dft_field_container() - for monitor_idx, monitor in enumerate(self._dft_convergence_monitors): - for component_idx, component in enumerate( - self.dft_field_components): - previous_fields = self._previous_fields[monitor_idx][ - component_idx] - current_fields[monitor_idx][component_idx] = sim.get_dft_array( - monitor, - component, - int(monitor.nfreqs // 2), - ) - norm_previous = _norm_fn(previous_fields) - field_diff = previous_fields - current_fields[monitor_idx][ - component_idx] - if norm_previous != 0: - relative_change.append( - _norm_fn(field_diff) / norm_previous) - else: - relative_change.append(1.0) - relative_change = _reduce_fn(relative_change) - self._dft_relative_change.append(relative_change) - self._previous_fields = current_fields - if mp.am_master() and mp.verbosity > 0: - print( - 'At simulation time %.2f the relative change in the DFT fields is %.2e.' - % (sim.meep_time(), relative_change)) - return relative_change < self.dft_threshold - - def _simulation_run_callback(self, sim: mp.Simulation) -> bool: - """A callback function returning `True` when the simulation should stop. - - This is a step function that gets called at each time step of the Meep - simulation, taking a Meep simulation object as an input and returning `True` - when the simulation should be terminated and returning `False` when the - simulation should continue. The resulting step function is designed to be - used with the `until` and `until_after_sources` arguments of the Meep - simulation `run()` routine. - - Args: - sim: a Meep simulation object. - - Returns: - a boolean indicating whether the simulation should be terminated. - """ - current_meep_time = sim.round_time() - if current_meep_time <= self._last_measurement_meep_time + self.measurement_interval: - return False - if current_meep_time <= self.minimum_run_time: - if mp.am_master() and mp.verbosity > 0: - remaining_time = self.minimum_run_time - sim.round_time() - self._log_fn( - '%.2f to go until the minimum simulation runtime is reached.' - % (remaining_time, )) - return False - if current_meep_time >= self.maximum_run_time: - if mp.am_master() and mp.verbosity > 0: - self._log_fn( - 'Stopping the simulation because the maximum simulation run ' - 'time of %.2f has been reached.' % - (self.maximum_run_time, )) - return True - self._last_measurement_meep_time = current_meep_time - return self._are_dfts_converged(sim) diff --git a/python/meep.i b/python/meep.i index 59ac569d7..ca47a0062 100644 --- a/python/meep.i +++ b/python/meep.i @@ -844,7 +844,9 @@ meep::volume_list *make_volume_list(const meep::volume &v, int c, //-------------------------------------------------- %inline %{ -void _get_gradient(PyObject *grad, PyObject *fields_a, PyObject *fields_f, PyObject *grid_volume, PyObject *frequencies, PyObject *py_geom_list, PyObject *f, bool sim_is_cylindrical) { +void _get_gradient(PyObject *grad, double scalegrad, PyObject *fields_a, PyObject *fields_f, + meep::grid_volume *grid_volume, meep::volume *where, PyObject *frequencies, + meep_geom::geom_epsilon *geps, PyObject *fields_shapes) { // clean the gradient array PyArrayObject *pao_grad = (PyArrayObject *)grad; if (!PyArray_Check(pao_grad)) meep::abort("grad parameter must be numpy array."); @@ -867,15 +869,11 @@ void _get_gradient(PyObject *grad, PyObject *fields_a, PyObject *fields_f, PyObj if (PyArray_NDIM(pao_fields_f) !=1) {meep::abort("Numpy forward fields array must have 1 dimension.");} std::complex *fields_f_c = (std::complex *)PyArray_DATA(pao_fields_f); - // scalegrad not currently used - double scalegrad = 1.0; - - // clean the volume object - void* where; - - PyObject *swigobj = PyObject_GetAttrString(grid_volume, "swigobj"); - SWIG_ConvertPtr(swigobj,&where,0,NULL); - const meep::volume* where_vol = (const meep::volume*)where; + // clean shapes array + PyArrayObject *pao_fields_shapes = (PyArrayObject *)fields_shapes; + if (!PyArray_Check(pao_fields_shapes)) meep::abort("fields shape parameter must be numpy array."); + if (!PyArray_ISCARRAY(pao_fields_shapes)) meep::abort("Numpy fields shape array must be C-style contiguous."); + size_t *fields_shapes_c = (size_t *)PyArray_DATA(pao_fields_shapes); // clean the frequencies array PyArrayObject *pao_freqs = (PyArrayObject *)frequencies; @@ -885,24 +883,8 @@ void _get_gradient(PyObject *grad, PyObject *fields_a, PyObject *fields_f, PyObj 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); - // prepare a geometry_tree - // TODO eventually it would be nice to store the geom tree within the structure object so we don't have to recreate it here - geometric_object_list *l; - l = new geometric_object_list(); - if (!py_list_to_gobj_list(py_geom_list,l)) meep::abort("Unable to convert geometry tree."); - geom_box_tree geometry_tree = calculate_tree(*where_vol,*l); - - // clean the fields pointer - void* f_v; - SWIG_ConvertPtr(f,&f_v,NULL,NULL); - meep::fields* f_c = (meep::fields *)f_v; - // calculate the gradient - meep_geom::material_grids_addgradient(grad_c,ng,fields_a_c,fields_f_c,frequencies_c,nf,scalegrad,*where_vol,geometry_tree,f_c, sim_is_cylindrical); - - destroy_geom_box_tree(geometry_tree); - delete l; - Py_DECREF(swigobj); + meep_geom::material_grids_addgradient(grad_c,ng,fields_a_c,fields_f_c,fields_shapes_c,frequencies_c,scalegrad,*grid_volume,*where,geps); } %} @@ -1413,11 +1395,6 @@ void _get_gradient(PyObject *grad, PyObject *fields_a, PyObject *fields_f, PyObj // For some reason SWIG needs the namespaced version too %apply material_type_list { meep_geom::material_type_list }; -// Typemaps for geom_epsilon object -%typemap(out) geom_epsilon { - $result = PyCapsule_New($1); -} - // Typemap suite for kpoint_func %typecheck(SWIG_TYPECHECK_POINTER) (meep::kpoint_func user_kpoint_func, void *user_kpoint_data) { diff --git a/python/tests/test_adjoint_solver.py b/python/tests/test_adjoint_solver.py index e65a6a842..84ea42768 100644 --- a/python/tests/test_adjoint_solver.py +++ b/python/tests/test_adjoint_solver.py @@ -15,6 +15,7 @@ resolution = 30 silicon = mp.Medium(epsilon=12) +sapphire = mp.Medium(epsilon_diag=(10.225,10.225,9.95),epsilon_offdiag=(-0.825,-0.55*np.sqrt(3/2),0.55*np.sqrt(3/2))) sxy = 5.0 cell_size = mp.Vector3(sxy,sxy,0) @@ -60,10 +61,10 @@ k_point = mp.Vector3(0.23,-0.38) -def forward_simulation(design_params, mon_type, frequencies=None): +def forward_simulation(design_params, mon_type, frequencies=None, mat2=silicon): matgrid = mp.MaterialGrid(mp.Vector3(Nx,Ny), mp.air, - silicon, + mat2, weights=design_params.reshape(Nx,Ny)) matgrid_geometry = [mp.Block(center=mp.Vector3(), @@ -114,10 +115,10 @@ def forward_simulation(design_params, mon_type, frequencies=None): return Ez2 -def adjoint_solver(design_params, mon_type, frequencies=None): +def adjoint_solver(design_params, mon_type, frequencies=None, mat2=silicon): matgrid = mp.MaterialGrid(mp.Vector3(Nx,Ny), mp.air, - silicon, + mat2, weights=np.ones((Nx,Ny))) matgrid_region = mpa.DesignRegion(matgrid, @@ -399,6 +400,31 @@ def test_complex_fields(self): tol = 0.005 if mp.is_single_precision() else 0.0008 self.assertClose(adj_scale,fd_grad,epsilon=tol) + def test_offdiagonal(self): + print("*** TESTING OFFDIAGONAL COMPONENTS ***") + + ## test the single frequency and multi frequency case + for frequencies in [[fcen], [1/1.58, fcen, 1/1.53]]: + ## compute gradient using adjoint solver + adjsol_obj, adjsol_grad = adjoint_solver(p, MonitorObject.EIGENMODE, frequencies, sapphire) + ## compute unperturbed S12 + S12_unperturbed = forward_simulation(p, MonitorObject.EIGENMODE, frequencies, sapphire) + + ## compare objective results + print("S12 -- adjoint solver: {}, traditional simulation: {}".format(adjsol_obj,S12_unperturbed)) + self.assertClose(adjsol_obj,S12_unperturbed,epsilon=1e-6) + + ## compute perturbed S12 + S12_perturbed = forward_simulation(p+dp, MonitorObject.EIGENMODE, frequencies, sapphire) + + ## compare gradients + if adjsol_grad.ndim < 2: + adjsol_grad = np.expand_dims(adjsol_grad,axis=1) + adj_scale = (dp[None,:]@adjsol_grad).flatten() + fd_grad = S12_perturbed-S12_unperturbed + print("Directional derivative -- adjoint solver: {}, FD: {}".format(adj_scale,fd_grad)) + tol = 0.04 + self.assertClose(adj_scale,fd_grad,epsilon=tol) if __name__ == '__main__': unittest.main() diff --git a/src/loop_in_chunks.cpp b/src/loop_in_chunks.cpp index eabd4c825..3c9bba764 100644 --- a/src/loop_in_chunks.cpp +++ b/src/loop_in_chunks.cpp @@ -229,7 +229,7 @@ void chunkloop_field_components::update_values(ptrdiff_t idx) { component of equal_shift (which should be either -2, 0, or +2). (equal_shift is there to prevent us from counting edge points twice.) */ -static ivec vec2diel_floor(const vec &pt, double a, const ivec &equal_shift) { +ivec fields::vec2diel_floor(const vec &pt, double a, const ivec &equal_shift) { ivec ipt(pt.dim); LOOP_OVER_DIRECTIONS(pt.dim, d) { ipt.set_direction(d, 1 + 2 * int(floor(pt.in_direction(d) * a - .5))); @@ -238,7 +238,7 @@ static ivec vec2diel_floor(const vec &pt, double a, const ivec &equal_shift) { } return ipt; } -static ivec vec2diel_ceil(const vec &pt, double a, const ivec &equal_shift) { +ivec fields::vec2diel_ceil(const vec &pt, double a, const ivec &equal_shift) { ivec ipt(pt.dim); LOOP_OVER_DIRECTIONS(pt.dim, d) { ipt.set_direction(d, 1 + 2 * int(ceil(pt.in_direction(d) * a - .5))); diff --git a/src/material_data.cpp b/src/material_data.cpp index 323c72eb4..6bf043a23 100644 --- a/src/material_data.cpp +++ b/src/material_data.cpp @@ -87,6 +87,7 @@ void material_data::copy_from(const material_data &from) { beta = from.beta; eta = from.eta; damping = from.damping; + material_grid_kinds = from.material_grid_kinds; } material_type_list::material_type_list() : items(NULL), num_items(0) {} diff --git a/src/material_data.hpp b/src/material_data.hpp index 2fcbae81d..94cea0d4f 100644 --- a/src/material_data.hpp +++ b/src/material_data.hpp @@ -136,6 +136,7 @@ struct material_data { double beta; double eta; double damping; + bool trivial=true; /* There are several possible scenarios when material grids overlap -- these different scenarios enable different applications. diff --git a/src/meep.hpp b/src/meep.hpp index ad77fb5e4..488c517b9 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -1931,6 +1931,8 @@ class fields { int is_phasing(); // loop_in_chunks.cpp + static ivec vec2diel_floor(const vec &pt, double a, const ivec &equal_shift); + static ivec vec2diel_ceil(const vec &pt, double a, const ivec &equal_shift); void loop_in_chunks(field_chunkloop chunkloop, void *chunkloop_data, const volume &where, component cgrid = Centered, bool use_symmetry = true, bool snap_unit_dims = false); diff --git a/src/meep/vec.hpp b/src/meep/vec.hpp index ec9b0f81b..34ca5cc98 100644 --- a/src/meep/vec.hpp +++ b/src/meep/vec.hpp @@ -789,6 +789,12 @@ class ivec { return result; }; + ivec operator+(int s) const { + ivec result = *this; + LOOP_OVER_DIRECTIONS(dim, d) result.t[d] += s; + return result; + }; + ivec operator+=(const ivec &a) { LOOP_OVER_DIRECTIONS(dim, d) t[d] += a.t[d]; return *this; @@ -847,6 +853,19 @@ class ivec { return result; }; + // element-wise product + ivec operator*(const ivec &a) const { + ivec result = *this; + LOOP_OVER_DIRECTIONS(dim, d) result.t[d] *= a.t[d]; + return result; + }; + + ivec operator/(int s) const { + ivec result = *this; + LOOP_OVER_DIRECTIONS(dim, d) result.t[d] /= s; + return result; + }; + vec operator*(double s) const { vec result(dim); LOOP_OVER_DIRECTIONS(dim, d) result.set_direction(d, t[d] * s); diff --git a/src/meepgeom.cpp b/src/meepgeom.cpp index c914e8097..1f57f052b 100644 --- a/src/meepgeom.cpp +++ b/src/meepgeom.cpp @@ -562,7 +562,7 @@ void epsilon_material_grid(material_data *md, double u) { // Linearly interpolate dc epsilon values cinterp_tensors(m1->epsilon_diag, m1->epsilon_offdiag, m2->epsilon_diag, m2->epsilon_offdiag, &mm->epsilon_diag, &mm->epsilon_offdiag, u); - + // Interpolate resonant strength from d.p. vector3 zero_vec; zero_vec.x = zero_vec.y = zero_vec.z = 0; @@ -604,11 +604,16 @@ void epsilon_material_grid(material_data *md, double u) { // TODO: dampen the lorentzians to improve stability // mm->D_conductivity_diag.x = mm->D_conductivity_diag.y = mm->D_conductivity_diag.z = u*(1-u) * // omega_mean; + md->trivial=false; } double fake_damping = u*(1-u)*(md->damping); mm->D_conductivity_diag.x += fake_damping; mm->D_conductivity_diag.y += fake_damping; mm->D_conductivity_diag.z += fake_damping; + + // set the trivial flag + if (md->damping != 0) + md->trivial=false; } // return material of the point p from the file (assumed already read) @@ -636,7 +641,16 @@ void epsilon_file_material(material_data *md, vector3 p) { geom_epsilon::geom_epsilon(geometric_object_list g, material_type_list mlist, const meep::volume &v) { - geometry = g; // don't bother making a copy, only used in one place + // Copy the geometry + int length = g.num_items; + geometry.num_items = length; + geometry.items = new geometric_object[length]; + for (int i = 0; i < length; i++){ + geometric_object_copy(&g.items[i],&geometry.items[i]); + geometry.items[i].material = new material_data(); + static_cast(geometry.items[i].material)->copy_from(*(material_data *)(g.items[i].material)); + } + extra_materials = mlist; current_pol = NULL; @@ -681,7 +695,16 @@ geom_epsilon::geom_epsilon(geometric_object_list g, material_type_list mlist, // copy constructor geom_epsilon::geom_epsilon(const geom_epsilon &geps1) { - geometry = geps1.geometry; // TODO make a copy + // Copy the geometry + int length = geps1.geometry.num_items; + geometry.num_items = length; + geometry.items = new geometric_object[length]; + for (int i = 0; i < length; i++){ + geometric_object_copy(&geps1.geometry.items[i],&geometry.items[i]); + geometry.items[i].material = new material_data(); + static_cast(geometry.items[i].material)->copy_from(*(material_data *)(geps1.geometry.items[i].material)); + } + geometry_tree = geps1.geometry_tree; restricted_tree = geps1.restricted_tree; extra_materials = geps1.extra_materials; @@ -691,6 +714,11 @@ geom_epsilon::geom_epsilon(const geom_epsilon &geps1) { } geom_epsilon::~geom_epsilon() { + int length = geometry.num_items; + for (int i = 0; i < length; i++){ + delete geometry.items[i].material; + } + delete[] geometry.items; unset_volume(); destroy_geom_box_tree(geometry_tree); FOR_DIRECTIONS(d) FOR_SIDES(b) { @@ -815,7 +843,7 @@ void geom_epsilon::get_material_pt(material_type &material, const meep::vec &r) tp = geom_tree_search(p, restricted_tree, &oi); // interpolate and project onto material grid - u = tanh_projection(matgrid_val(p, tp, oi, md), md->beta, md->eta); + u = tanh_projection(matgrid_val(p, tp, oi, md)+this->u_p, md->beta, md->eta); // interpolate material from material grid point epsilon_material_grid(md, u); @@ -1279,7 +1307,7 @@ void geom_epsilon::fallback_chi1inv_row(meep::component c, double chi1inv_row[3] int oi; tp = geom_tree_search(p, restricted_tree, &oi); gradient = matgrid_grad(p, tp, oi, md); - uval = matgrid_val(p, tp, oi, md); + uval = matgrid_val(p, tp, oi, md)+this->u_p; } else { gradient = normal_vector(meep::type(c), v); @@ -1622,7 +1650,7 @@ void geom_epsilon::sigma_row(meep::component c, double sigrow[3], const meep::ve tp = geom_tree_search(p, restricted_tree, &oi); // interpolate and project onto material grid - u = tanh_projection(matgrid_val(p, tp, oi, mat), mat->beta, mat->eta); + u = tanh_projection(matgrid_val(p, tp, oi, mat)+this->u_p, mat->beta, mat->eta); epsilon_material_grid(mat, u); // interpolate material from material grid point mat->medium.check_offdiag_im_zero_or_abort(); } @@ -1963,6 +1991,11 @@ set the materials as specified */ void set_materials_from_geom_epsilon(meep::structure *s, geom_epsilon *geps, vector3 center, bool use_anisotropic_averaging, double tol, int maxeval, absorber_list alist) { + + // store for later use in gradient calculations + geps->tol = tol; + geps->maxeval = maxeval; + meep::grid_volume gv = s->gv; if (alist) { for (absorber_list_type::iterator layer = alist->begin(); layer != alist->end(); layer++) { @@ -2503,20 +2536,43 @@ double vec_to_value(vector3 diag, vector3 offdiag, int idx) { return val; } -void get_material_tensor(const medium_struct *mm, double freq, - std::complex *tensor) { - /* - Note that the current implementation assumes that any dispersion - is either a lorentzian or drude profile. - */ +void invert_tensor(std::complex t_inv[9], std::complex t[9]) { + +#define m(x,y) t[x*3+y] +#define minv(x,y) t_inv[x*3+y] + std::complex det = m(0, 0) * (m(1, 1) * m(2, 2) - m(2, 1) * m(1, 2)) - + m(0, 1) * (m(1, 0) * m(2, 2) - m(1, 2) * m(2, 0)) + + m(0, 2) * (m(1, 0) * m(2, 1) - m(1, 1) * m(2, 0)); + std::complex invdet = 1.0 / det; + minv(0, 0) = (m(1, 1) * m(2, 2) - m(2, 1) * m(1, 2)) * invdet; + minv(0, 1) = (m(0, 2) * m(2, 1) - m(0, 1) * m(2, 2)) * invdet; + minv(0, 2) = (m(0, 1) * m(1, 2) - m(0, 2) * m(1, 1)) * invdet; + minv(1, 0) = (m(1, 2) * m(2, 0) - m(1, 0) * m(2, 2)) * invdet; + minv(1, 1) = (m(0, 0) * m(2, 2) - m(0, 2) * m(2, 0)) * invdet; + minv(1, 2) = (m(1, 0) * m(0, 2) - m(0, 0) * m(1, 2)) * invdet; + minv(2, 0) = (m(1, 0) * m(2, 1) - m(2, 0) * m(1, 1)) * invdet; + minv(2, 1) = (m(2, 0) * m(0, 1) - m(0, 0) * m(2, 1)) * invdet; + minv(2, 2) = (m(0, 0) * m(1, 1) - m(1, 0) * m(0, 1)) * invdet; +#undef m(x,y) +#undef minv(x,y) +} + +void eff_chi1inv_row_disp(meep::component c, std::complex chi1inv_row[3], + const meep::vec &r, double freq, geom_epsilon *geps) { + // locate the proper material + material_type md; + geps->get_material_pt(md, r); + const medium_struct *mm = &(md->medium); + std::complex tensor[9], tensor_inv[9]; + + // loop over all the tensor components for (int i = 0; i < 9; i++) { - std::complex a = std::complex(1, 0); - std::complex b = std::complex(0, 0); + std::complex a,b; // compute first part containing conductivity vector3 dummy; dummy.x = dummy.y = dummy.z = 0.0; double conductivityCur = vec_to_value(mm->D_conductivity_diag, dummy, i); - a += std::complex(0.0, conductivityCur / freq); + a = std::complex(1.0, conductivityCur / (freq)); // compute lorentzian component b = cvec_to_value(mm->epsilon_diag, mm->epsilon_offdiag, i); @@ -2530,66 +2586,104 @@ void get_material_tensor(const medium_struct *mm, double freq, // elementwise multiply tensor[i] = a * b; } -} - -double get_material_gradient( - double u, // material parameter at current point - std::complex fields_a, // adjoint field at current point - std::complex fields_f, // forward field at current point - double freq, // frequency - material_data *md, // material - meep::component field_dir, // current field component - double du // step size -) { - /* Note that the current implementation assumes that - no materials have off-diag components and that if a material - has dispersion, it is either a lorentzian or drude profile. - */ - const medium_struct *mm = &(md->medium); - const medium_struct *m1 = &(md->medium_1); - const medium_struct *m2 = &(md->medium_2); - - // trivial case - if ((mm->E_susceptibilities.size() == 0) && mm->D_conductivity_diag.x == 0 && - mm->D_conductivity_diag.y == 0 && mm->D_conductivity_diag.z == 0){ - switch (field_dir){ - case meep::Er: - case meep::Ex: return (m2->epsilon_diag.x - m1->epsilon_diag.x) * (fields_a * fields_f).real(); - case meep::Ep: - case meep::Ey: return (m2->epsilon_diag.y - m1->epsilon_diag.y) * (fields_a * fields_f).real(); - case meep::Ez: return (m2->epsilon_diag.z - m1->epsilon_diag.z) * (fields_a * fields_f).real(); - default: meep::abort("Invalid field component specified in gradient."); - } - } + // invert the matrix + invert_tensor(tensor_inv, tensor); - /* For now we do a finite difference approach to estimate the - gradient of the system matrix `A` since it's fairly accurate, - cheap, and easy to generalize. */ - std::complex dA_du_0[9] = {std::complex(0, 0)}; - epsilon_material_grid(md, u - du); - get_material_tensor(mm, freq, dA_du_0); - - std::complex dA_du_1[9] = {std::complex(0, 0)}; - epsilon_material_grid(md, u + du); - get_material_tensor(mm, freq, dA_du_1); + // get the row we care about + switch (component_direction(c)) { + case meep::X: + case meep::R: + chi1inv_row[0] = tensor_inv[0]; + chi1inv_row[1] = tensor_inv[1]; + chi1inv_row[2] = tensor_inv[2]; + break; + case meep::Y: + case meep::P: + chi1inv_row[0] = tensor_inv[3]; + chi1inv_row[1] = tensor_inv[4]; + chi1inv_row[2] = tensor_inv[5]; + break; + case meep::Z: + chi1inv_row[0] = tensor_inv[6]; + chi1inv_row[1] = tensor_inv[7]; + chi1inv_row[2] = tensor_inv[8]; + break; + case meep::NO_DIRECTION: chi1inv_row[0] = chi1inv_row[1] = chi1inv_row[2] = 0; break; + } +} - std::complex dA_du[9] = {std::complex(0, 0)}; - for (int i = 0; i < 9; i++) - dA_du[i] = (dA_du_1[i] - dA_du_0[i]) / (2 * du); +std::complex get_material_gradient( + const meep::vec &r, // current point + const meep::component adjoint_c, // adjoint field component + const meep::component forward_c, // forward field component + std::complex fields_f, // forward field at current point + double freq, // frequency + geom_epsilon *geps, // material + meep::grid_volume &gv, // simulation grid volume + double du // step size +) { + /*Compute the Aᵤx product from the -λᵀAᵤx calculation. + The current adjoint (λ) field component (adjoint_c) + determines which row of Aᵤ we care about. + The current forward (x) field component (forward_c) + determines which column of Aᵤ we care about. + + There are two ways we can compute the required row/ + column of Aᵤ: + 1. If we want to incorporate subpixel smoothing, + we use the eff_chi1inv_row() routine. This neglects + conductivities, susceptibilities, etc. + 2. We use eff_chi1inv_row_disp() for all other cases + (at the expense of not accounting for subpixel smoothing, + if there were any). + + For now we do a finite difference approach to estimate the + gradient of the system matrix A since it's fairly accurate, + cheap, and easy to generalize.*/ + + // locate the proper material + material_type md; + geps->get_material_pt(md, r); int dir_idx = 0; - if (field_dir == meep::Ex || field_dir == meep::Er) + if (forward_c == meep::Ex || forward_c == meep::Er) dir_idx = 0; - else if (field_dir == meep::Ey || field_dir == meep::Ep) + else if (forward_c == meep::Ey || forward_c == meep::Ep) dir_idx = 1; - else if (field_dir == meep::Ez) + else if (forward_c == meep::Ez) dir_idx = 2; else meep::abort("Invalid adjoint field component"); - std::complex result = fields_a * dA_du[3 * dir_idx + dir_idx] * fields_f; - return result.real(); + if (md->trivial){ + const double sd = 1.0; // FIXME: make user-changable? + meep::volume v(r); + LOOP_OVER_DIRECTIONS(dim, d){ + v.set_direction_min(d, r.in_direction(d) - 0.5*gv.inva*sd); + v.set_direction_max(d, r.in_direction(d) + 0.5*gv.inva*sd); + } + double row_1[3], row_2[3], dA_du[3]; + geps->u_p = -du; + geps->eff_chi1inv_row(adjoint_c, row_1, v, geps->tol, geps->maxeval); + geps->u_p = du; + geps->eff_chi1inv_row(adjoint_c, row_2, v, geps->tol, geps->maxeval); + geps->u_p=0; + + for (int i=0;i<3;i++) dA_du[i] = (row_1[i] - row_2[i])/(2*du); + return dA_du[dir_idx] * fields_f; + + }else{ + std::complex row_1[3], row_2[3], dA_du[3]; + geps->u_p = -du; + eff_chi1inv_row_disp(adjoint_c,row_1,r,freq,geps); + geps->u_p = du; + eff_chi1inv_row_disp(adjoint_c,row_2,r,freq,geps); + geps->u_p=0; + + for (int i=0;i<3;i++) dA_du[i] = (row_1[i] - row_2[i])/(2*du); + return dA_du[dir_idx] * fields_f; + } } /* add the weights from linear_interpolate (see the linear_interpolate @@ -2635,16 +2729,15 @@ in row-major order (the order used by HDF5): */ #undef U } -void material_grids_addgradient_point(double *v, std::complex fields_a, - std::complex fields_f, meep::component field_dir, - vector3 p, double scalegrad, double freq, - geom_box_tree geometry_tree) { +void material_grids_addgradient_point(double *v, + vector3 p, double scalegrad, + geom_epsilon *geps) { geom_box_tree tp; int oi, ois; material_data *mg, *mg_sum; double uval; int kind; - tp = geom_tree_search(p, geometry_tree, &oi); + tp = geom_tree_search(p, geps->geometry_tree, &oi); if (tp && ((material_type)tp->objects[oi].o->material)->which_subclass == material_data::MATERIAL_GRID) @@ -2658,7 +2751,7 @@ void material_grids_addgradient_point(double *v, std::complex fields_a, if ((tp) && ((kind = mg->material_grid_kinds) == material_data::U_MEAN)) { int matgrid_val_count = 0; geom_box_tree tp_sum; - tp_sum = geom_tree_search(p, geometry_tree, &ois); + tp_sum = geom_tree_search(p, geps->geometry_tree, &ois); mg_sum = (material_data *)tp_sum->objects[ois].o->material; do { tp_sum = geom_tree_search_next(p, tp_sum, &ois); @@ -2688,12 +2781,12 @@ void material_grids_addgradient_point(double *v, std::complex fields_a, vector3 sz = mg->grid_size; double *vcur = v; double *ucur = mg->weights; - uval = matgrid_val(p, tp, oi, mg); + uval = tanh_projection(matgrid_val(p, tp, oi, mg), mg->beta, mg->eta); 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, - get_material_gradient(uval, fields_a, fields_f, freq, mg, field_dir, 1e-6) * scalegrad, + scalegrad, ucur, kind, uval); if (kind == material_data::U_DEFAULT) break; tp = geom_tree_search_next(p, tp, &oi); @@ -2707,99 +2800,181 @@ void material_grids_addgradient_point(double *v, std::complex fields_a, 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, - get_material_gradient(uval, fields_a, fields_f, freq, mg, field_dir) * - scalegrad, + scalegrad, ucur, kind, uval); } } -void material_grids_addgradient(double *v, size_t ng, std::complex *fields_a, - std::complex *fields_f, double *frequencies, - size_t nf, double scalegrad, const meep::volume &where, - geom_box_tree geometry_tree, meep::fields *f, bool sim_is_cylindrical) { - size_t n2, n3, n4; - double s[3][3], cen[3][3], c1, c2, c3, s1, s2, s3; - vector3 p; +/**********************************************************/ +/* Some gradient helper routines */ +/**********************************************************/ - // calculate cell dimensions - meep::direction dirs[3]; - meep::vec min_max_loc[2] = {meep::vec(0,0,0),meep::vec(0,0,0)}; // extremal points in subgrid - meep::component field_dir[3] = {meep::Ex, meep::Ey, meep::Ez}; - if (sim_is_cylindrical) { - field_dir[0] = meep::Er; - field_dir[1] = meep::Ep; - } - size_t dims[9] = {1, 1, 1, 1, 1, 1, 1, 1, 1}; - for (int c = 0; c < 3; c++) { - - f->get_array_slice_dimensions(where, &dims[3 * c], dirs, true, false, min_max_loc, - 0, field_dir[c]); - - vector3 max_corner = vec_to_vector3(min_max_loc[1]); - double max_c_array[3] = {max_corner.x, max_corner.y, max_corner.z}; - vector3 min_corner = vec_to_vector3(min_max_loc[0]); - double min_c_array[3] = {min_corner.x, min_corner.y, min_corner.z}; - - if (sim_is_cylindrical){ - max_c_array[0] = max_corner.z; - max_c_array[1] = max_corner.x; - max_c_array[2] = max_corner.y; - min_c_array[0] = min_corner.z; - min_c_array[1] = min_corner.x; - min_c_array[2] = min_corner.y; - } +#define LOOP_OVER_DIRECTIONS_BACKWARDS(dim, d) \ + for (meep::direction d = (meep::direction)(meep::stop_at_direction(dim)-1), \ + loop_stop_directi = meep::start_at_direction(dim); \ + d >= loop_stop_directi; d = (meep::direction)(d - 1)) - for (int ci = 0; ci < 3; ci++) { - s[c][ci] = (max_c_array[ci] - min_c_array[ci]) == 0 ? 0 : (max_c_array[ci] - min_c_array[ci]) / (dims[3 * c + ci] - 1); - cen[c][ci] = dims[3 * c + ci] <= 1 ? 0 : min_c_array[ci]; +void set_strides(meep::ndim dim, ptrdiff_t *the_stride,const meep::ivec c1,const meep::ivec c2) { + FOR_DIRECTIONS(d) { the_stride[d] = 1;} + meep::ivec n_s = (c2 - c1)/2 + 1; + LOOP_OVER_DIRECTIONS_BACKWARDS(dim,d) { + ptrdiff_t current_stride = 1; + LOOP_OVER_DIRECTIONS_BACKWARDS(dim,d_i) { + if (d_i==d){ + the_stride[d] = current_stride; + break; + } + current_stride *= n_s.in_direction(d_i); } + } +} + +ptrdiff_t get_idx_from_ivec(meep::ndim dim, meep::ivec c1, ptrdiff_t *the_stride, meep::ivec v){ + ptrdiff_t idx = 0; + meep::ivec diff = ((v-c1) / 2); + LOOP_OVER_DIRECTIONS(dim,d){ + idx += diff.in_direction(d)*the_stride[d]; } + return idx; +} - // Loop over component, x, y, z, and frequency dimensions - // TODO speed up with MPI (if needed) - int xyz_index = 0; - if (sim_is_cylindrical){ - for (int c = 0; c < 3; c++) { // component - n2 = dims[c * 3]; n3 = dims[c * 3 + 1]; n4 = dims[c * 3 + 2]; - c1 = cen[c][0]; c2 = cen[c][1]; c3 = cen[c][2]; - s1 = s[c][0]; s2 = s[c][1]; s3 = s[c][2]; - - for (size_t i1 = 0; i1 < nf; ++i1) { // freq - for (size_t i2 = 0; i2 < n2; ++i2) { // z - for (size_t i4 = 0; i4 < n4; ++i4) {//y - for (size_t i3 = 0; i3 < n3; ++i3) {//x - p.z = i2 * s1 + c1; p.x = i3 * s2 + c2; p.y = i4 * s3 + c3; - material_grids_addgradient_point(v+ ng*i1, fields_a[xyz_index]*p.x, fields_f[xyz_index], field_dir[c], p, - scalegrad, frequencies[i1], geometry_tree); - //p.x is the (2pi)r' factor from integrating in cyldrical coordinate; - //2pi is canceled out by a overcouted factor of 2pi*r of the Near2FarRegion; See near2far.cpp - xyz_index++; - } - } - } - } +void material_grids_addgradient(double *v, size_t ng, std::complex *fields_a, + std::complex *fields_f, size_t fields_shapes[12], double *frequencies, + double scalegrad, meep::grid_volume &gv, + meep::volume &where, geom_epsilon *geps) { + +// only loop over field components relevant to our simulation (in the proper order) +#define FOR_MY_COMPONENTS(c1) FOR_ELECTRIC_COMPONENTS(c1) if (!coordinate_mismatch(gv.dim, component_direction(c1))) + + /* poach some logic from loop_in_chunks + that ensures we loop over the same grid + points that the DFT points lie on. */ + meep::ivec *is = new meep::ivec[3]; + meep::ivec *ie = new meep::ivec[3]; + int c_i = 0; + FOR_MY_COMPONENTS(cgrid) { + meep::vec yee_c(gv.yee_shift(meep::Centered) - gv.yee_shift(cgrid)); + meep::ivec iyee_c(gv.iyee_shift(meep::Centered) - gv.iyee_shift(cgrid)); + meep::volume wherec(where + yee_c); + is[c_i] = meep::ivec(meep::fields::vec2diel_floor(wherec.get_min_corner(), gv.a, zero_ivec(gv.dim))); + ie[c_i] = meep::ivec(meep::fields::vec2diel_ceil(wherec.get_max_corner(), gv.a, zero_ivec(gv.dim))); + c_i++; + } + //calculate the number of elements in an entire block (x,y,z) for each component + size_t nf = fields_shapes[0]; + size_t stride_row[3] = {1,1,1}; //first entry is a dummy entry to simplify indexing + for (int i=0;i<3;i++) { + for (int j=1;j<4;j++){ + stride_row[i] *= fields_shapes[4*i+j]; } - } else { - for (int c = 0; c < 3; c++) { // component - n2 = dims[c * 3]; n3 = dims[c * 3 + 1]; n4 = dims[c * 3 + 2]; - c1 = cen[c][0]; c2 = cen[c][1]; c3 = cen[c][2]; - s1 = s[c][0]; s2 = s[c][1]; s3 = s[c][2]; - - for (size_t i1 = 0; i1 < nf; ++i1) { // freq - for (size_t i2 = 0; i2 < n2; ++i2) { // x - for (size_t i3 = 0; i3 < n3; ++i3) { // y - for (size_t i4 = 0; i4 < n4; ++i4) { // z - p.x = i2 * s1 + c1; p.y = i3 * s2 + c2; p.z = i4 * s3 + c3; - material_grids_addgradient_point(v+ ng*i1, fields_a[xyz_index], fields_f[xyz_index], field_dir[c], p, - scalegrad, frequencies[i1], geometry_tree); - xyz_index++; + } + size_t c_start[3]; + c_start[0] = 0; + c_start[1] = nf*stride_row[0]; + c_start[2] = nf*(stride_row[0]+stride_row[1]); + +// fields dimensions are (components, nfreqs, x, y, z) +#define GET_FIELDS(fields,comp,freq,idx) fields[c_start[comp] + freq*stride_row[comp] + idx] + + meep::ivec start_ivec; + // loop over frequency + for (size_t f_i = 0; f_i < nf; ++f_i) { + int ci_adjoint = 0; + FOR_MY_COMPONENTS(adjoint_c) { + LOOP_OVER_IVECS(gv, is[ci_adjoint], ie[ci_adjoint], idx) { + size_t idx_fields = IVEC_LOOP_COUNTER; + meep::ivec ip = gv.iloc(adjoint_c,idx); + meep::vec p = gv.loc(adjoint_c,idx); + std::complex adj = GET_FIELDS(fields_a,ci_adjoint,f_i,idx_fields); + double cyl_scale; + int ci_forward = 0; + FOR_MY_COMPONENTS(forward_c) { + /* we need to calculate the bounds of + the forward fields (in space) so that we + can properly index into the fields array + later */ + meep::ivec isf = is[ci_forward]; + meep::ivec ief = ie[ci_forward]; + // the actual starting index, as if we were running another LOOP_OVER_IVECS + ptrdiff_t idx0_f = (isf - (gv).little_corner()).yucky_val(0) / 2 * loop_s1 + \ + (isf - (gv).little_corner()).yucky_val(1) / 2 * loop_s2 + \ + (isf - (gv).little_corner()).yucky_val(2) / 2 * loop_s3; + start_ivec = gv.iloc(forward_c,idx0_f); + ptrdiff_t the_stride[5]; + set_strides(gv.dim, the_stride, isf,ief); + /**************************************/ + /* Main Routine */ + /**************************************/ + /********* compute -λᵀAᵤx *************/ + + /* trivial case, no interpolation/restriction needed */ + if (forward_c == adjoint_c){ + std::complex fwd = GET_FIELDS(fields_f,ci_forward,f_i,idx_fields); + std::complex prod = adj * get_material_gradient(p, adjoint_c, forward_c, fwd, frequencies[f_i], geps, gv, 1e-6); + cyl_scale = (gv.dim == meep::Dcyl) ? 2*p.r() : 1; // the pi is already factored in near2far.cpp + material_grids_addgradient_point(v+ng*f_i, vec_to_vector3(p), scalegrad*prod.real()*cyl_scale, geps); + /* more complicated case requires interpolation/restriction */ + }else{ + /* we need to restrict the adjoint fields to + the two nodes of interest, and interpolate + the forward fields to the same two nodes. Then + we perform our inner product at these nodes + */ + std::complex fwd_avg, fwd1, fwd2, prod; + ptrdiff_t fwd1_idx, fwd2_idx; + + //identify the first corner of the forward fields + meep::ivec fwd_p = ip + gv.iyee_shift(forward_c) - gv.iyee_shift(adjoint_c); + + //identify the other three corners + 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 fwd_pa = (fwd_p + unit_a*2); + meep::ivec fwd_pf = (fwd_p - unit_f*2); + meep::ivec fwd_paf = (fwd_p + unit_a*2 - unit_f*2); + + //identify the two eps points + meep::ivec ieps1 = (fwd_p + fwd_pf) / 2; + meep::ivec ieps2 = (fwd_pa + fwd_paf) / 2; + + //operate on the first eps node + fwd1_idx = get_idx_from_ivec(gv.dim, start_ivec, the_stride, fwd_p); + fwd1 = 0.5 * GET_FIELDS(fields_f,ci_forward,f_i,fwd1_idx); + fwd2_idx = get_idx_from_ivec(gv.dim, start_ivec, the_stride, fwd_pf); + fwd2 = 0.5 * GET_FIELDS(fields_f,ci_forward,f_i,fwd2_idx); + fwd_avg = fwd1 + fwd2; + meep::vec eps1 = gv[ieps1]; + cyl_scale = (gv.dim == meep::Dcyl) ? eps1.r() : 1; + prod = 0.5*adj*get_material_gradient(eps1, adjoint_c, forward_c, fwd_avg, frequencies[f_i], geps, gv, 1e-6); + material_grids_addgradient_point(v+ng*f_i, vec_to_vector3(eps1), scalegrad*prod.real()*cyl_scale, geps); + + //operate on the second eps node + fwd1_idx = get_idx_from_ivec(gv.dim, start_ivec, the_stride, fwd_pa); + fwd1 = 0.5 * GET_FIELDS(fields_f,ci_forward,f_i,fwd1_idx); + fwd2_idx = get_idx_from_ivec(gv.dim, start_ivec, the_stride, fwd_paf); + fwd2 = 0.5 * GET_FIELDS(fields_f,ci_forward,f_i,fwd2_idx); + fwd_avg = fwd1 + fwd2; + meep::vec eps2 = gv[ieps2]; + cyl_scale = (gv.dim == meep::Dcyl) ? eps2.r() : 1; + prod = 0.5*adj*get_material_gradient(eps2, adjoint_c, forward_c, fwd_avg, frequencies[f_i], geps, gv, 1e-6); + material_grids_addgradient_point(v+ng*f_i, vec_to_vector3(eps2), scalegrad*prod.real()*cyl_scale, geps); } + ci_forward++; } + /********* compute λᵀbᵤ ***************/ + /* not yet implemented/needed */ + /**************************************/ + /**************************************/ } + ci_adjoint++; } } +#undef GET_FIELDS +#undef FOR_MY_COMPONENTS + delete[] is; + delete[] ie; } -} + static void find_array_min_max(int n, const double *data, double &min_val, double &max_val) { min_val = data[0]; diff --git a/src/meepgeom.hpp b/src/meepgeom.hpp index bd723866e..8995e7c63 100644 --- a/src/meepgeom.hpp +++ b/src/meepgeom.hpp @@ -170,13 +170,16 @@ struct cond_profile { }; class geom_epsilon : public meep::material_function { - geometric_object_list geometry; - geom_box_tree geometry_tree; - geom_box_tree restricted_tree; - - cond_profile cond[5][2]; // [direction][side] public: + double u_p = 0; + geom_box_tree geometry_tree; + geom_box_tree restricted_tree; + geometric_object_list geometry; + cond_profile cond[5][2]; // [direction][side] + double tol=DEFAULT_SUBPIXEL_TOL; + int maxeval=DEFAULT_SUBPIXEL_MAXEVAL; + geom_epsilon(geometric_object_list g, material_type_list mlist, const meep::volume &v); geom_epsilon(const geom_epsilon &geps1); // copy constructor virtual ~geom_epsilon(); @@ -214,9 +217,8 @@ class geom_epsilon : public meep::material_function { void add_susceptibilities(meep::structure *s); void add_susceptibilities(meep::field_type ft, meep::structure *s); -private: void get_material_pt(material_type &material, const meep::vec &r); - +private: material_type_list extra_materials; pol *current_pol; }; @@ -293,22 +295,10 @@ 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 get_material_tensor(const medium_struct *mm, double freq, std::complex *tensor); -double get_material_gradient(double u, std::complex fields_a, - std::complex fields_f, double freq, - material_data *md, meep::component field_dir, - double du = 1.0e-3); -void add_interpolate_weights(double rx, double ry, double rz, - double *data, int nx, int ny, int nz, int stride, - double scaleby, const double *udata, int ukind, double uval); -void material_grids_addgradient_point(double *v, std::complex fields_a, - std::complex fields_f, meep::component field_dir, - vector3 p, double scalegrad, double freq, - geom_box_tree geometry_tree); void material_grids_addgradient(double *v, size_t ng, std::complex *fields_a, - std::complex *fields_f, double *frequencies, - size_t nf, double scalegrad, const meep::volume &where, - geom_box_tree geometry_tree, meep::fields *f, bool sim_is_cylindrical); + std::complex *fields_f, size_t fields_shapes[12], + double *frequencies, double scalegrad, + meep::grid_volume &gv, meep::volume &where, geom_epsilon *geps); /***************************************************************/ /* routines in GDSIIgeom.cc ************************************/