diff --git a/doc/docs/Python_Tutorials/Adjoint_Solver.md b/doc/docs/Python_Tutorials/Adjoint_Solver.md index 1de1b7103..665817bca 100644 --- a/doc/docs/Python_Tutorials/Adjoint_Solver.md +++ b/doc/docs/Python_Tutorials/Adjoint_Solver.md @@ -204,286 +204,315 @@ runtime of this script using three Intel Xeon 2.0 GHz processors is approximately 14 hours. ```py -import numpy as np -import matplotlib -matplotlib.use("agg") -import matplotlib.pyplot as plt +from typing import List, NamedTuple, Tuple + from autograd import numpy as npa, tensor_jacobian_product, grad -import nlopt +import matplotlib.pyplot as plt import meep as mp import meep.adjoint as mpa -from typing import NamedTuple - -resolution = 50 # pixels/μm - -w = 0.4 # waveguide width -l = 3.0 # waveguide length (on each side of design region) -dpad = 0.6 # padding length above/below design region -dpml = 1.0 # PML thickness -dx = 1.6 # length of design region -dy = 1.6 # width of design region - -sx = dpml + l + dx + l + dpml -sy = dpml + dpad + dy + dpad + dpml -cell_size = mp.Vector3(sx, sy, 0) - -pml_layers = [mp.PML(thickness=dpml)] - -# wavelengths for minimax optimization -wvls = (1.265, 1.270, 1.275, 1.285, 1.290, 1.295) -frqs = [1 / wvl for wvl in wvls] - -minimum_length = 0.15 # minimum length scale (μm) -eta_i = 0.5 # blueprint design field thresholding point (between 0 and 1) -eta_e = 0.75 # erosion design field thresholding point (between 0 and 1) -eta_d = 1 - eta_e # dilation design field thresholding point (between 0 and 1) -filter_radius = mpa.get_conic_radius_from_eta_e(minimum_length, eta_e) -print(f"filter_radius:, {filter_radius:.6f}") - -# pulsed source center frequency and bandwidth -wvl_min = 1.26 -wvl_max = 1.30 -frq_min = 1 / wvl_max -frq_max = 1 / wvl_min -fcen = 0.5 * (frq_min + frq_max) -df = frq_max - frq_min - -eig_parity = mp.ODD_Z -src_pt = mp.Vector3(-0.5 * sx + dpml, 0, 0) - -nSiO2 = 1.5 -SiO2 = mp.Medium(index=nSiO2) -nSi = 3.5 -Si = mp.Medium(index=nSi) - -design_region_size = mp.Vector3(dx, dy, 0) -design_region_resolution = int(2 * resolution) -Nx = int(design_region_size.x * design_region_resolution) + 1 -Ny = int(design_region_size.y * design_region_resolution) + 1 - -# impose a bit "mask" of thickness equal to the filter radius -# around the edges of the design region in order to prevent -# violations of the minimum feature size constraint. - -x_g = np.linspace( - -design_region_size.x / 2, - design_region_size.x / 2, - Nx, -) -y_g = np.linspace( - -design_region_size.y / 2, - design_region_size.y / 2, - Ny, -) -X_g, Y_g = np.meshgrid( - x_g, - y_g, - sparse=True, - indexing="ij", -) +import nlopt +import numpy as np -left_wg_mask = (X_g <= -design_region_size.x / 2 + filter_radius) & ( - np.abs(Y_g) <= w / 2 + +RESOLUTION_UM = 50 +WAVELENGTH_MIN_UM = 1.26 +WAVELENGTH_MAX_UM = 1.30 +WAVEGUIDE_UM = mp.Vector3(0.4, 3.0, 0) +PADDING_UM = 0.6 +PML_UM = 1.0 +DESIGN_WAVELENGTHS_UM = (1.265, 1.270, 1.275, 1.285, 1.290, 1.295) +DESIGN_REGION_UM = mp.Vector3(1.6, 1.6, mp.inf) +DESIGN_REGION_RESOLUTION_UM = int(2 * RESOLUTION_UM) +NX_DESIGN_GRID = int(DESIGN_REGION_UM.x * DESIGN_REGION_RESOLUTION_UM) + 1 +NY_DESIGN_GRID = int(DESIGN_REGION_UM.y * DESIGN_REGION_RESOLUTION_UM) + 1 +MIN_LENGTH_UM = 0.15 +SIGMOID_THRESHOLD_INTRINSIC = 0.5 +SIGMOID_THRESHOLD_EROSION = 0.75 +SIGMOID_THRESHOLD_DILATION = 1 - SIGMOID_THRESHOLD_EROSION +MODE_SYMMETRY = mp.ODD_Z +SILICON = mp.Medium(index=3.5) +SILICON_DIOXIDE = mp.Medium(index=1.5) + +cell_um = mp.Vector3( + PML_UM + WAVEGUIDE_UM.x + DESIGN_REGION_UM.x + WAVEGUIDE_UM.x + PML_UM, + PML_UM + PADDING_UM + DESIGN_REGION_UM.y + PADDING_UM + PML_UM, + 0 ) -right_wg_mask = (X_g >= design_region_size.x / 2 - filter_radius) & ( - np.abs(Y_g) <= w / 2 +filter_radius_um = mpa.get_conic_radius_from_eta_e( + MIN_LENGTH_UM, + SIGMOID_THRESHOLD_EROSION ) -Si_mask = left_wg_mask | right_wg_mask +frequency_min = 1 / WAVELENGTH_MAX_UM +frequency_max = 1 / WAVELENGTH_MIN_UM +frequency_center = 0.5 * (frequency_min + frequency_max) +frequency_width = frequency_max - frequency_min +pml_layers = [mp.PML(thickness=PML_UM)] +frequencies = [1 / wavelength_um for wavelength_um in DESIGN_WAVELENGTHS_UM] +src_pt = mp.Vector3(-0.5 * cell_um.x + PML_UM, 0, 0) +refl_pt = mp.Vector3(-0.5 * cell_um.x + PML_UM + 0.5 * WAVEGUIDE_UM.x) +tran_pt = mp.Vector3(0.5 * cell_um.x - PML_UM - 0.5 * WAVEGUIDE_UM.x) +stop_cond = mp.stop_when_fields_decayed(50, mp.Ez, refl_pt, 1e-6) -border_mask = ( - (X_g <= -design_region_size.x / 2 + filter_radius) - | (X_g >= design_region_size.x / 2 - filter_radius) - | (Y_g <= -design_region_size.y / 2 + filter_radius) - | (Y_g >= design_region_size.y / 2 - filter_radius) -) -SiO2_mask = border_mask.copy() -SiO2_mask[Si_mask] = False -refl_pt = mp.Vector3(-0.5 * sx + dpml + 0.5 * l) -tran_pt = mp.Vector3(0.5 * sx - dpml - 0.5 * l) +def str_from_list(list_: List[float]) -> str: + return "[" + ", ".join(f"{val:.4f}" for val in list_) + "]" -stop_cond = mp.stop_when_fields_decayed(50, mp.Ez, refl_pt, 1e-8) +def border_masks() -> Tuple[np.ndarray, np.ndarray]: + """Return border masks for the design region. -def mapping(x: np.ndarray, eta: float, beta: float) -> np.ndarray: - """A differentiable mapping function which applies, in order, - the following sequence of transformations to the design weights: - (1) a bit mask for the boundary pixels, (2) convolution with a - conic filter, and (3) projection via a hyperbolic tangent (if - necessary). + The masks are used to prevent violations on constraints on the + minimum feature size at the boundaries of the design region. + + Returns: + A 2-tuple of 2D arrays for border masks for Si and SiO2. + """ + x_grid = np.linspace( + -DESIGN_REGION_UM.x / 2, + DESIGN_REGION_UM.x / 2, + NX_DESIGN_GRID, + ) + y_grid = np.linspace( + -DESIGN_REGION_UM.y / 2, + DESIGN_REGION_UM.y / 2, + NY_DESIGN_GRID, + ) + xy_grid_x, xy_grid_y = np.meshgrid( + x_grid, + y_grid, + sparse=True, + indexing="ij", + ) + + left_waveguide_port = ( + (xy_grid_x <= -DESIGN_REGION_UM.x / 2 + filter_radius_um) & + (np.abs(xy_grid_y) <= WAVEGUIDE_UM.y / 2) + ) + right_waveguide_port = ( + (xy_grid_x >= DESIGN_REGION_UM.x / 2 - filter_radius_um) & + (np.abs(xy_grid_y) <= WAVEGUIDE_UM.y / 2) + ) + silicon_mask = left_waveguide_port | right_waveguide_port + + border_mask = ( + (xy_grid_x <= -DESIGN_REGION_UM.x / 2 + filter_radius_um) | + (xy_grid_x >= DESIGN_REGION_UM.x / 2 - filter_radius_um) | + (xy_grid_y <= -DESIGN_REGION_UM.y / 2 + filter_radius_um) | + (xy_grid_y >= DESIGN_REGION_UM.y / 2 - filter_radius_um) + ) + silicon_dioxide_mask = border_mask.copy() + silicon_dioxide_mask[silicon_mask] = False + + return silicon_mask, silicon_dioxide_mask + + +def filter_projection( + weights: np.ndarray, + sigmoid_threshold: float, + sigmoid_bias: float +) -> np.ndarray: + """A differentiable function to filter and project the design weights. Args: - x: design weights as a 1d array of size Nx*Ny. - eta: erosion/dilation parameter for the projection. - beta: bias parameter for the projection. A value of 0 is no projection. + weights: design weights as a flattened (1D) array. + sigmoid_threshold: erosion/dilation parameter for the projection. + sigmoid_bias: bias parameter for the projection. 0 is no projection. Returns: - The mapped design weights as a 1d array. + The mapped design weights as a 1D array. """ - x = npa.where( - Si_mask.flatten(), + silicon_mask, silicon_dioxide_mask = border_masks() + + weights_masked = npa.where( + silicon_mask.flatten(), 1, npa.where( - SiO2_mask.flatten(), + silicon_dioxide_mask.flatten(), 0, - x, + weights, ), ) - filtered_field = mpa.conic_filter( - x, - filter_radius, - design_region_size.x, - design_region_size.y, - design_region_resolution, + weights_filtered = mpa.conic_filter( + weights_masked, + filter_radius_um, + DESIGN_REGION_UM.x, + DESIGN_REGION_UM.y, + DESIGN_REGION_RESOLUTION_UM, ) - if beta == 0: - return filtered_field.flatten() - + if sigmoid_bias == 0: + return weights_filtered.flatten() else: - projected_field = mpa.tanh_projection( - filtered_field, - beta, - eta, + weights_projected = mpa.tanh_projection( + weights_filtered, + sigmoid_bias, + sigmoid_threshold, ) - - return projected_field.flatten() + return weights_projected.flatten() -def f(x: np.ndarray, grad: np.ndarray) -> float: +def obj_func(epigraph_and_weights: np.ndarray, grad: np.ndarray) -> float: """Objective function for the epigraph formulation. Args: - x: 1d array of size 1+Nx*Ny containing epigraph variable (first element) - and design weights (remaining Nx*Ny elements). - grad: the gradient as a 1d array of size 1+Nx*Ny modified in place. + epigraph_and_weights: 1D array containing epigraph variable (first + element) and design weights (remaining elements). + grad: the gradient as a flattened (1D) array, modified in place. Returns: - The epigraph variable (a scalar). + The scalar epigraph variable. """ - t = x[0] # epigraph variable - v = x[1:] # design weights + epigraph = epigraph_and_weights[0] + if grad.size > 0: grad[0] = 1 grad[1:] = 0 - return t + + return epigraph -def c(result: np.ndarray, x: np.ndarray, gradient: np.ndarray, eta: float, - beta: float, use_epsavg: bool): +def epigraph_constraint( + result: np.ndarray, + epigraph_and_weights: np.ndarray, + gradient: np.ndarray, + sigmoid_threshold: float, + sigmoid_bias: float, + use_epsavg: bool, +) -> None: """Constraint function for the epigraph formulation. Args: - result: the result of the function evaluation modified in place. - x: 1d array of size 1+Nx*Ny containing epigraph variable (first - element) and design weights (remaining Nx*Ny elements). - gradient: the Jacobian matrix with dimensions (1+Nx*Ny, - 2*num. wavelengths) modified in place. - eta: erosion/dilation parameter for projection. - beta: bias parameter for projection. + result: the result of the function evaluation, modified in place. + epigraph_and_weights: 1D array containing the epigraph variable (first + element) and design weights (remaining elements). + gradient: the Jacobian matrix with dimensions (1 + NX_DESIGN_GRID * + NY_DESIGN_GRID, 2 * num. wavelengths), modified in place. + sigmoid_threshold: erosion/dilation parameter for projection. + sigmoid_bias: bias parameter for projection. use_epsavg: whether to use subpixel smoothing. """ - t = x[0] # epigraph variable - v = x[1:] # design weights - - f0, dJ_du = opt([mapping(v, eta, 0 if use_epsavg else beta)]) + epigraph = epigraph_and_weights[0] + weights = epigraph_and_weights[1:] + + obj_val, grad = opt( + [ + filter_projection( + weights, + sigmoid_threshold, + 0 if use_epsavg else sigmoid_bias + ) + ] + ) - f0_reflection = f0[0] - f0_transmission = f0[1] - f0_merged = np.concatenate((f0_reflection, f0_transmission)) - f0_merged_str = '[' + ','.join(str(ff) for ff in f0_merged) + ']' + reflectance = obj_val[0] + transmittance = obj_val[1] + obj_val_merged = np.concatenate((reflectance, transmittance)) + obj_val_merged_str = str_from_list(obj_val_merged) - dJ_du_reflection = dJ_du[0] - dJ_du_transmission = dJ_du[1] - nfrq = len(frqs) - my_grad = np.zeros((Nx * Ny, 2 * nfrq)) - my_grad[:, :nfrq] = dJ_du_reflection - my_grad[:, nfrq:] = dJ_du_transmission + grad_reflectance = grad[0] + grad_transmittance = grad[1] + nfrq = len(frequencies) + grad = np.zeros((NX_DESIGN_GRID * NY_DESIGN_GRID, 2 * nfrq)) + grad[:, :nfrq] = grad_reflectance + grad[:, nfrq:] = grad_transmittance - # backpropagate the gradients through mapping function + # backpropagate the gradients through filter and projection function for k in range(2 * nfrq): - my_grad[:, k] = tensor_jacobian_product(mapping, 0)( - v, - eta, - beta, - my_grad[:, k], + grad[:, k] = tensor_jacobian_product(filter_projection, 0)( + weights, + sigmoid_threshold, + sigmoid_bias, + grad[:, k], ) if gradient.size > 0: - gradient[:, 0] = -1 # gradient w.r.t. epigraph variable ("t") - gradient[:, 1:] = my_grad.T # gradient w.r.t. each frequency objective + gradient[:, 0] = -1 # gradient w.r.t. epigraph variable + gradient[:, 1:] = grad.T # gradient w.r.t. each frequency objective - result[:] = np.real(f0_merged) - t + result[:] = np.real(obj_val_merged) - epigraph - objfunc_history.append(np.real(f0_merged)) - epivar_history.append(t) + objfunc_history.append(np.real(obj_val_merged)) + epivar_history.append(epigraph) print( - f"iteration:, {cur_iter[0]:3d}, eta: {eta}, beta: {beta:2d}, " - f"t: {t:.5f}, obj. func.: {f0_merged_str}" + f"iteration:, {cur_iter[0]:3d}, sigmoid_bias: {sigmoid_bias:2d}, " + f"epigraph: {epigraph:.5f}, obj. func.: {obj_val_merged_str}" ) cur_iter[0] = cur_iter[0] + 1 -def glc(result: np.ndarray, x: np.ndarray, gradient: np.ndarray, - beta: float) -> float: +def line_width_and_spacing_constraint( + result: np.ndarray, + epigraph_and_weights: np.ndarray, + gradient: np.ndarray, + sigmoid_bias: float +) -> float: """Constraint function for the minimum linewidth. Args: result: the result of the function evaluation modified in place. - x: 1d array of size 1+Nx*Ny containing epigraph variable (first - element) and design weights (remaining elements). - gradient: the Jacobian matrix with dimensions (1+Nx*Ny, - num. wavelengths) modified in place. - beta: bias parameter for projection. + epigraph_and_weights: 1D array containing epigraph variable (first + element) and design weights (remaining elements). + gradient: the Jacobian matrix, modified in place. + sigmoid_bias: bias parameter for projection. Returns: The value of the constraint function (a scalar). """ - t = x[0] # dummy parameter - v = x[1:] # design parameters + epigraph = epigraph_and_weights[0] + weights = epigraph_and_weights[1:] a1 = 1e-3 # hyper parameter (primary) b1 = 0 # hyper parameter (secondary) gradient[:, 0] = -a1 - filter_f = lambda a: mpa.conic_filter( - a.reshape(Nx, Ny), - filter_radius, - design_region_size.x, - design_region_size.y, - design_region_resolution, + filter_func = lambda a: mpa.conic_filter( + a.reshape(NX_DESIGN_GRID, NY_DESIGN_GRID), + filter_radius_um, + DESIGN_REGION_UM.x, + DESIGN_REGION_UM.y, + DESIGN_REGION_RESOLUTION_UM, + ) + + threshold_func = lambda a: mpa.tanh_projection( + a, + sigmoid_bias, + SIGMOID_THRESHOLD_INTRINSIC ) - threshold_f = lambda a: mpa.tanh_projection(a, beta, eta_i) # hyper parameter (constant factor and exponent) - c0 = 1e7 * (filter_radius * 1 / resolution) ** 4 + c0 = 1e7 * (filter_radius_um * 1 / RESOLUTION_UM) ** 4 - M1 = lambda a: mpa.constraint_solid(a, c0, eta_e, filter_f, threshold_f, 1) - M2 = lambda a: mpa.constraint_void(a, c0, eta_d, filter_f, threshold_f, 1) + M1 = lambda a: mpa.constraint_solid( + a, c0, SIGMOID_THRESHOLD_EROSION, filter_func, threshold_func, 1 + ) + M2 = lambda a: mpa.constraint_void( + a, c0, SIGMOID_THRESHOLD_DILATION, filter_func, threshold_func, 1 + ) - g1 = grad(M1)(v) - g2 = grad(M2)(v) + g1 = grad(M1)(weights) + g2 = grad(M2)(weights) - result[0] = M1(v) - a1 * t - b1 - result[1] = M2(v) - a1 * t - b1 + result[0] = M1(weights) - a1 * epigraph - b1 + result[1] = M2(weights) - a1 * epigraph - b1 gradient[0, 1:] = g1.flatten() gradient[1, 1:] = g2.flatten() - t1 = (M1(v) - b1) / a1 - t2 = (M2(v) - b1) / a1 + t1 = (M1(weights) - b1) / a1 + t2 = (M2(weights) - b1) / a1 - print(f"glc:, {result[0]}, {result[1]}, {t1}, {t2}") + print(f"line_width_and_spacing_constraint:, {result[0]}, {result[1]}, {t1}, {t2}") return max(t1, t2) def straight_waveguide() -> (np.ndarray, NamedTuple): - """Computes the DFT fields from the mode source in a straight waveguide - for use as normalization of the reflectance measurement during the - optimization. + """Computes the DFT fields from the mode source in a straight waveguide. + + The DFT fields are used as normalization of the reflectance measurement + during the optimization. Returns: A 2-tuple consisting of a 1d array of DFT fields and DFT fields object @@ -491,26 +520,26 @@ def straight_waveguide() -> (np.ndarray, NamedTuple): """ sources = [ mp.EigenModeSource( - src=mp.GaussianSource(fcen, fwidth=df), - size=mp.Vector3(0, sy, 0), + src=mp.GaussianSource(frequency_center, fwidth=frequency_width), + size=mp.Vector3(0, cell_um.y, 0), center=src_pt, eig_band=1, - eig_parity=eig_parity, + eig_parity=MODE_SYMMETRY, ) ] geometry = [ mp.Block( - size=mp.Vector3(mp.inf, w, mp.inf), + size=mp.Vector3(mp.inf, WAVEGUIDE_UM.y, mp.inf), center=mp.Vector3(), - material=Si, + material=SILICON, ) ] sim = mp.Simulation( - resolution=resolution, - default_material=SiO2, - cell_size=cell_size, + resolution=RESOLUTION_UM, + default_material=SILICON_DIOXIDE, + cell_size=cell_um, sources=sources, geometry=geometry, boundary_layers=pml_layers, @@ -518,8 +547,8 @@ def straight_waveguide() -> (np.ndarray, NamedTuple): ) refl_mon = sim.add_mode_monitor( - frqs, - mp.ModeRegion(center=refl_pt, size=mp.Vector3(0, sy, 0)), + frequencies, + mp.ModeRegion(center=refl_pt, size=mp.Vector3(0, cell_um.y, 0)), yee_grid=True, ) @@ -528,7 +557,7 @@ def straight_waveguide() -> (np.ndarray, NamedTuple): res = sim.get_eigenmode_coefficients( refl_mon, [1], - eig_parity=eig_parity, + eig_parity=MODE_SYMMETRY, ) coeffs = res.alpha @@ -539,15 +568,16 @@ def straight_waveguide() -> (np.ndarray, NamedTuple): def mode_converter_optimization( - input_flux: np.ndarray, - input_flux_data: NamedTuple, - use_damping: bool, - use_epsavg: bool, - beta: float) -> mpa.OptimizationProblem: + input_flux: np.ndarray, + input_flux_data: NamedTuple, + use_damping: bool, + use_epsavg: bool, + sigmoid_bias: float, +) -> mpa.OptimizationProblem: """Sets up the adjoint optimization of the waveguide mode converter. Args: - input_flux: 1d array of DFT fields from normalization run. + input_flux: 1D array of DFT fields from the normalization run. input_flux_data: DFT fields object returned by `meep.get_flux_data`. use_damping: whether to use the damping feature of `MaterialGrid`. use_epsavg: whether to use subpixel smoothing in `MaterialGrid`. @@ -556,20 +586,20 @@ def mode_converter_optimization( A `meep.adjoint.OptimizationProblem` class object. """ matgrid = mp.MaterialGrid( - mp.Vector3(Nx, Ny, 0), - SiO2, - Si, - weights=np.ones((Nx, Ny)), - beta=beta if use_epsavg else 0, + mp.Vector3(NX_DESIGN_GRID, NY_DESIGN_GRID, 0), + SILICON_DIOXIDE, + SILICON, + weights=np.ones((NX_DESIGN_GRID, NY_DESIGN_GRID)), + beta=sigmoid_bias if use_epsavg else 0, do_averaging=True if use_epsavg else False, - damping=0.02 * 2 * np.pi * fcen if use_damping else 0, + damping=0.02 * 2 * np.pi * frequency_center if use_damping else 0, ) matgrid_region = mpa.DesignRegion( matgrid, volume=mp.Volume( center=mp.Vector3(), - size=mp.Vector3(design_region_size.x, design_region_size.y, mp.inf), + size=DESIGN_REGION_UM ), ) @@ -584,8 +614,8 @@ def mode_converter_optimization( geometry = [ mp.Block( center=mp.Vector3(), - size=mp.Vector3(mp.inf, w, mp.inf), - material=Si, + size=mp.Vector3(mp.inf, WAVEGUIDE_UM.y, mp.inf), + material=SILICON, ) ] @@ -593,18 +623,18 @@ def mode_converter_optimization( sources = [ mp.EigenModeSource( - src=mp.GaussianSource(fcen, fwidth=df), - size=mp.Vector3(0, sy, 0), + src=mp.GaussianSource(frequency_center, fwidth=frequency_width), + size=mp.Vector3(0, cell_um.y, 0), center=src_pt, eig_band=1, - eig_parity=eig_parity, + eig_parity=MODE_SYMMETRY, ), ] sim = mp.Simulation( - resolution=resolution, - default_material=SiO2, - cell_size=cell_size, + resolution=RESOLUTION_UM, + default_material=SILICON_DIOXIDE, + cell_size=cell_um, sources=sources, geometry=geometry, boundary_layers=pml_layers, @@ -616,21 +646,21 @@ def mode_converter_optimization( sim, mp.Volume( center=refl_pt, - size=mp.Vector3(0, sy, 0), + size=mp.Vector3(0, cell_um.y, 0), ), 1, forward=False, - eig_parity=eig_parity, + eig_parity=MODE_SYMMETRY, subtracted_dft_fields=input_flux_data, ), mpa.EigenmodeCoefficient( sim, mp.Volume( center=tran_pt, - size=mp.Vector3(0, sy, 0), + size=mp.Vector3(0, cell_um.y, 0), ), 2, - eig_parity=eig_parity, + eig_parity=MODE_SYMMETRY, ), ] @@ -645,7 +675,7 @@ def mode_converter_optimization( objective_functions=[J1, J2], objective_arguments=obj_list, design_regions=[matgrid_region], - frequencies=frqs, + frequencies=frequencies, ) return opt @@ -654,114 +684,136 @@ def mode_converter_optimization( if __name__ == "__main__": input_flux, input_flux_data = straight_waveguide() - algorithm = nlopt.LD_MMA + silicon_mask, silicon_dioxide_mask = border_masks() - # number of design parameters - n = Nx * Ny + num_weights = NX_DESIGN_GRID * NY_DESIGN_GRID - # initial design parameters - x = np.ones((n,)) * 0.5 - x[Si_mask.flatten()] = 1.0 # set the edges of waveguides to silicon - x[SiO2_mask.flatten()] = 0.0 # set the other edges to SiO2 + # Initial design weights. + epigraph_and_weights = np.ones((num_weights,)) * 0.5 + epigraph_and_weights[silicon_mask.flatten()] = 1.0 + epigraph_and_weights[silicon_dioxide_mask.flatten()] = 0.0 - # lower and upper bounds for design weights - lb = np.zeros((n,)) - lb[Si_mask.flatten()] = 1.0 - ub = np.ones((n,)) - ub[SiO2_mask.flatten()] = 0.0 + # Lower and upper bounds for design weights. + weights_lower_bound = np.zeros((num_weights,)) + weights_lower_bound[silicon_mask.flatten()] = 1.0 + weights_upper_bound = np.ones((num_weights,)) + weights_upper_bound[silicon_dioxide_mask.flatten()] = 0.0 - # insert epigraph variable initial value (arbitrary) and bounds into the - # design array. the actual value is determined by the objective and + # Insert epigraph variable with arbitrary initial value and bounds into the + # design array. The actual value is determined by the objective and # constraint functions below. - x = np.insert(x, 0, 1.2) - lb = np.insert(lb, 0, -np.inf) - ub = np.insert(ub, 0, +np.inf) + epigraph_and_weights = np.insert(epigraph_and_weights, 0, 1.2) + weights_lower_bound = np.insert(weights_lower_bound, 0, -np.inf) + weights_upper_bound = np.insert(weights_upper_bound, 0, +np.inf) objfunc_history = [] epivar_history = [] cur_iter = [0] - beta_thresh = 64 # threshold beta above which to use subpixel smoothing - betas = [8, 16, 32, 64, 128, 256] + sigmoid_bias_threshold = 64 # threshold beta above which to use subpixel smoothing + sigmoid_biases = [8, 16, 32, 64, 128, 256] max_evals = [80, 80, 100, 120, 120, 100] - tol_epi = np.array([1e-4] * 2 * len(frqs)) # R, 1-T - tol_lw = np.array([1e-8] * 2) # line width, line spacing - - for beta, max_eval in zip(betas, max_evals): - solver = nlopt.opt(algorithm, n + 1) - solver.set_lower_bounds(lb) - solver.set_upper_bounds(ub) - solver.set_min_objective(f) + tolerance_epigraph = np.array([1e-4] * 2 * len(frequencies)) # R, 1-T + tolerance_width_and_spacing = np.array([1e-8] * 2) # line width, line spacing + + for sigmoid_bias, max_eval in zip(sigmoid_biases, max_evals): + solver = nlopt.opt(nlopt.LD_MMA, num_weights + 1) + solver.set_lower_bounds(weights_lower_bound) + solver.set_upper_bounds(weights_upper_bound) + solver.set_min_objective(obj_func) solver.set_maxeval(max_eval) solver.set_param("dual_ftol_rel", 1e-7) solver.add_inequality_mconstraint( - lambda rr, xx, gg: c( - rr, - xx, - gg, - eta_i, - beta, - False if beta < beta_thresh else True, + lambda result_, epigraph_and_weights_, grad_: epigraph_constraint( + result_, + epigraph_and_weights_, + grad_, + SIGMOID_THRESHOLD_INTRINSIC, + sigmoid_bias, + False if sigmoid_bias < sigmoid_bias_threshold else True, ), - tol_epi, + tolerance_epigraph, ) solver.set_param("verbosity", 1) + if sigmoid_bias < sigmoid_bias_threshold: + use_epsavg = False + else: + use_epsavg = True + opt = mode_converter_optimization( input_flux, input_flux_data, True, # use_damping - False if beta < beta_thresh else True, # use_epsavg - beta, + use_epsavg, + sigmoid_bias, ) - # apply the minimum linewidth constraint + # Apply the minimum linewidth constraint # only in the final epoch to an initial # binary design from the previous epoch. - if beta == betas[-1]: - res = np.zeros(2) - grd = np.zeros((2, n + 1)) - t = glc(res, x, grd, beta) + if sigmoid_bias == sigmoid_biases[-1]: + line_width_and_spacing = np.zeros(2) + grad_line_width_and_spacing = np.zeros((2, num_weights + 1)) + linewidth_constraint_val = line_width_and_spacing_constraint( + line_width_and_spacing, + epigraph_and_weights, + grad_line_width_and_spacing, + sigmoid_bias + ) solver.add_inequality_mconstraint( - lambda rr, xx, gg: glc( - rr, - xx, - gg, - beta, + lambda result_, epigraph_and_weights_, grad_: line_width_and_spacing_constraint( + result_, + epigraph_and_weights_, + grad_, + sigmoid_bias, ), - tol_lw, + tolerance_width_and_spacing, ) - # execute a single forward run before the start of each + # Execute a single forward run before the start of each # epoch and manually set the initial epigraph variable to # slightly larger than the largest value of the objective # function over the six wavelengths and the lengthscale # constraint (final epoch only). - t0 = opt( + epigraph_initial = opt( [ - mapping( - x[1:], - eta_i, - beta if beta < beta_thresh else 0, + filter_projection( + epigraph_and_weights[1:], + SIGMOID_THRESHOLD_INTRINSIC, + sigmoid_bias if sigmoid_bias < sigmoid_bias_threshold else 0, ), ], need_gradient=False, ) - t0 = np.concatenate((t0[0][0], t0[0][1])) - t0_str = '[' + ','.join(str(tt) for tt in t0) + ']' - x[0] = np.amax(t0) - x[0] = 1.05 * (max(x[0], t) if beta == betas[-1] else x[0]) - print(f"data:, {beta}, {t0_str}, {x[0]}") + epigraph_initial = np.concatenate( + ( + epigraph_initial[0][0], + epigraph_initial[0][1] + ) + ) + epigraph_initial_str = str_from_list(epigraph_initial) - x[:] = solver.optimize(x) + epigraph_and_weights[0] = np.amax(epigraph_initial) + if sigmoid_bias == sigmoid_biases[-1]: + epigraph_and_weights[0] = 1.05 * max( + epigraph_and_weights[0], + linewidth_constraint_val + ) + print( + f"epigraph-calibration:, {sigmoid_bias}, {epigraph_initial_str}, " + f"{epigraph_and_weights[0]}" + ) + + epigraph_and_weights[:] = solver.optimize(epigraph_and_weights) - optimal_design_weights = mapping( - x[1:], - eta_i, - beta, - ).reshape(Nx, Ny) + optimal_design_weights = filter_projection( + epigraph_and_weights[1:], + SIGMOID_THRESHOLD_INTRINSIC, + sigmoid_bias, + ).reshape(NX_DESIGN_GRID, NY_DESIGN_GRID) - # save the unmapped weights and a bitmap image + # Save the unmapped weights and a bitmap image # of the design weights at the end of each epoch. fig, ax = plt.subplots() ax.imshow( @@ -772,34 +824,42 @@ if __name__ == "__main__": ax.set_axis_off() if mp.am_master(): fig.savefig( - f"optimal_design_beta{beta}.png", + f"optimal_design_beta{sigmoid_bias}.png", dpi=150, bbox_inches="tight", ) - # save the final (unmapped) design as a 2d array in CSV format + # Save the final (unmapped) design as a 2D array in CSV format np.savetxt( - f"unmapped_design_weights_beta{beta}.csv", - x[1:].reshape(Nx, Ny), + f"unmapped_design_weights_beta{sigmoid_bias}.csv", + epigraph_and_weights[1:].reshape( + NX_DESIGN_GRID, + NY_DESIGN_GRID + ), fmt="%4.2f", delimiter=",", ) - # save all the important optimization parameters and output - # as separate arrays in a single file for post processing. - with open("optimal_design.npz", "wb") as fl: + # Save important optimization parameters and output as + # separate fields for post processing. + with open("optimal_design.npz", "wb") as datafile: np.savez( - fl, - Nx=Nx, - Ny=Ny, - design_region_size=(dx, dy), - design_region_resolution=design_region_resolution, - betas=betas, + datafile, + RESOLUTION_UM=RESOLUTION_UM, + DESIGN_WAVELENGTHS_UM=DESIGN_WAVELENGTHS_UM, + WAVEGUIDE_UM=WAVEGUIDE_UM, + PADDING_UM=PADDING_UM, + PML_UM=PML_UM, + DESIGN_REGION_UM=DESIGN_REGION_UM, + DESIGN_REGION_RESOLUTION_UM=DESIGN_REGION_RESOLUTION_UM, + NX_DESIGN_GRID=NX_DESIGN_GRID, + NY_DESIGN_GRID=NY_DESIGN_GRID, + MIN_LENGTH_UM=MIN_LENGTH_UM, + sigmoid_biases=sigmoid_biases, max_eval=max_eval, objfunc_history=objfunc_history, epivar_history=epivar_history, - t=x[0], - unmapped_design_weights=x[1:], - minimum_length=minimum_length, + epigraph_variable=epigraph_and_weights[0], + unmapped_design_weights=epigraph_and_weights[1:], optimal_design_weights=optimal_design_weights, ) ``` diff --git a/python/examples/adjoint_optimization/mode_converter.py b/python/examples/adjoint_optimization/mode_converter.py index d33991949..b8543d4eb 100644 --- a/python/examples/adjoint_optimization/mode_converter.py +++ b/python/examples/adjoint_optimization/mode_converter.py @@ -1,6 +1,4 @@ -"""Topology optimization of the waveguide mode converter using - Meep's adjoint solver from A.M. Hammond et al., Optics Express, - Vol. 30, pp. 4467-4491 (2022). doi.org/10.1364/OE.442074 +"""Topology optimization of a waveguide mode converter. The worst-case optimization is based on minimizing the maximum of {R,1-T} where R (reflectance) is $|S_{11}|^2$ for mode 1 @@ -11,292 +9,306 @@ Vol. 29, pp. 23916-23938, (2021). doi.org/10.1364/OE.431188 """ -import numpy as np -import matplotlib +from typing import List, NamedTuple, Tuple -matplotlib.use("agg") -import matplotlib.pyplot as plt from autograd import numpy as npa, tensor_jacobian_product, grad -import nlopt +import matplotlib.pyplot as plt import meep as mp import meep.adjoint as mpa -from typing import NamedTuple - -resolution = 50 # pixels/μm - -w = 0.4 # waveguide width -l = 3.0 # waveguide length (on each side of design region) -dpad = 0.6 # padding length above/below design region -dpml = 1.0 # PML thickness -dx = 1.6 # length of design region -dy = 1.6 # width of design region - -sx = dpml + l + dx + l + dpml -sy = dpml + dpad + dy + dpad + dpml -cell_size = mp.Vector3(sx, sy, 0) - -pml_layers = [mp.PML(thickness=dpml)] - -# wavelengths for minimax optimization -wvls = (1.265, 1.270, 1.275, 1.285, 1.290, 1.295) -frqs = [1 / wvl for wvl in wvls] - -minimum_length = 0.15 # minimum length scale (μm) -eta_i = 0.5 # blueprint design field thresholding point (between 0 and 1) -eta_e = 0.75 # erosion design field thresholding point (between 0 and 1) -eta_d = 1 - eta_e # dilation design field thresholding point (between 0 and 1) -filter_radius = mpa.get_conic_radius_from_eta_e(minimum_length, eta_e) -print(f"filter_radius:, {filter_radius:.6f}") - -# pulsed source center frequency and bandwidth -wvl_min = 1.26 -wvl_max = 1.30 -frq_min = 1 / wvl_max -frq_max = 1 / wvl_min -fcen = 0.5 * (frq_min + frq_max) -df = frq_max - frq_min - -eig_parity = mp.ODD_Z -src_pt = mp.Vector3(-0.5 * sx + dpml, 0, 0) - -nSiO2 = 1.5 -SiO2 = mp.Medium(index=nSiO2) -nSi = 3.5 -Si = mp.Medium(index=nSi) - -design_region_size = mp.Vector3(dx, dy, 0) -design_region_resolution = int(2 * resolution) -Nx = int(design_region_size.x * design_region_resolution) + 1 -Ny = int(design_region_size.y * design_region_resolution) + 1 - -# impose a bit "mask" of thickness equal to the filter radius -# around the edges of the design region in order to prevent -# violations of the minimum feature size constraint. - -x_g = np.linspace( - -design_region_size.x / 2, - design_region_size.x / 2, - Nx, -) -y_g = np.linspace( - -design_region_size.y / 2, - design_region_size.y / 2, - Ny, -) -X_g, Y_g = np.meshgrid( - x_g, - y_g, - sparse=True, - indexing="ij", -) +import nlopt +import numpy as np -left_wg_mask = (X_g <= -design_region_size.x / 2 + filter_radius) & ( - np.abs(Y_g) <= w / 2 + +RESOLUTION_UM = 50 +WAVELENGTH_MIN_UM = 1.26 +WAVELENGTH_MAX_UM = 1.30 +WAVEGUIDE_UM = mp.Vector3(0.4, 3.0, 0) +PADDING_UM = 0.6 +PML_UM = 1.0 +DESIGN_WAVELENGTHS_UM = (1.265, 1.270, 1.275, 1.285, 1.290, 1.295) +DESIGN_REGION_UM = mp.Vector3(1.6, 1.6, mp.inf) +DESIGN_REGION_RESOLUTION_UM = int(2 * RESOLUTION_UM) +NX_DESIGN_GRID = int(DESIGN_REGION_UM.x * DESIGN_REGION_RESOLUTION_UM) + 1 +NY_DESIGN_GRID = int(DESIGN_REGION_UM.y * DESIGN_REGION_RESOLUTION_UM) + 1 +MIN_LENGTH_UM = 0.15 +SIGMOID_THRESHOLD_INTRINSIC = 0.5 +SIGMOID_THRESHOLD_EROSION = 0.75 +SIGMOID_THRESHOLD_DILATION = 1 - SIGMOID_THRESHOLD_EROSION +MODE_SYMMETRY = mp.ODD_Z +SILICON = mp.Medium(index=3.5) +SILICON_DIOXIDE = mp.Medium(index=1.5) + +cell_um = mp.Vector3( + PML_UM + WAVEGUIDE_UM.x + DESIGN_REGION_UM.x + WAVEGUIDE_UM.x + PML_UM, + PML_UM + PADDING_UM + DESIGN_REGION_UM.y + PADDING_UM + PML_UM, + 0, ) -right_wg_mask = (X_g >= design_region_size.x / 2 - filter_radius) & ( - np.abs(Y_g) <= w / 2 +filter_radius_um = mpa.get_conic_radius_from_eta_e( + MIN_LENGTH_UM, SIGMOID_THRESHOLD_EROSION ) -Si_mask = left_wg_mask | right_wg_mask +frequency_min = 1 / WAVELENGTH_MAX_UM +frequency_max = 1 / WAVELENGTH_MIN_UM +frequency_center = 0.5 * (frequency_min + frequency_max) +frequency_width = frequency_max - frequency_min +pml_layers = [mp.PML(thickness=PML_UM)] +frequencies = [1 / wavelength_um for wavelength_um in DESIGN_WAVELENGTHS_UM] +src_pt = mp.Vector3(-0.5 * cell_um.x + PML_UM, 0, 0) +refl_pt = mp.Vector3(-0.5 * cell_um.x + PML_UM + 0.5 * WAVEGUIDE_UM.x) +tran_pt = mp.Vector3(0.5 * cell_um.x - PML_UM - 0.5 * WAVEGUIDE_UM.x) +stop_cond = mp.stop_when_fields_decayed(50, mp.Ez, refl_pt, 1e-6) -border_mask = ( - (X_g <= -design_region_size.x / 2 + filter_radius) - | (X_g >= design_region_size.x / 2 - filter_radius) - | (Y_g <= -design_region_size.y / 2 + filter_radius) - | (Y_g >= design_region_size.y / 2 - filter_radius) -) -SiO2_mask = border_mask.copy() -SiO2_mask[Si_mask] = False -refl_pt = mp.Vector3(-0.5 * sx + dpml + 0.5 * l) -tran_pt = mp.Vector3(0.5 * sx - dpml - 0.5 * l) +def str_from_list(list_: List[float]) -> str: + return "[" + ", ".join(f"{val:.4f}" for val in list_) + "]" + + +def border_masks() -> Tuple[np.ndarray, np.ndarray]: + """Return border masks for the design region. + + The masks are used to prevent violations on constraints on the + minimum feature size at the boundaries of the design region. + + Returns: + A 2-tuple of 2D arrays for border masks for Si and SiO2. + """ + x_grid = np.linspace( + -DESIGN_REGION_UM.x / 2, + DESIGN_REGION_UM.x / 2, + NX_DESIGN_GRID, + ) + y_grid = np.linspace( + -DESIGN_REGION_UM.y / 2, + DESIGN_REGION_UM.y / 2, + NY_DESIGN_GRID, + ) + xy_grid_x, xy_grid_y = np.meshgrid( + x_grid, + y_grid, + sparse=True, + indexing="ij", + ) + + left_waveguide_port = (xy_grid_x <= -DESIGN_REGION_UM.x / 2 + filter_radius_um) & ( + np.abs(xy_grid_y) <= WAVEGUIDE_UM.y / 2 + ) + right_waveguide_port = (xy_grid_x >= DESIGN_REGION_UM.x / 2 - filter_radius_um) & ( + np.abs(xy_grid_y) <= WAVEGUIDE_UM.y / 2 + ) + silicon_mask = left_waveguide_port | right_waveguide_port + + border_mask = ( + (xy_grid_x <= -DESIGN_REGION_UM.x / 2 + filter_radius_um) + | (xy_grid_x >= DESIGN_REGION_UM.x / 2 - filter_radius_um) + | (xy_grid_y <= -DESIGN_REGION_UM.y / 2 + filter_radius_um) + | (xy_grid_y >= DESIGN_REGION_UM.y / 2 - filter_radius_um) + ) + silicon_dioxide_mask = border_mask.copy() + silicon_dioxide_mask[silicon_mask] = False -stop_cond = mp.stop_when_fields_decayed(50, mp.Ez, refl_pt, 1e-8) + return silicon_mask, silicon_dioxide_mask -def mapping(x: np.ndarray, eta: float, beta: float) -> np.ndarray: - """A differentiable mapping function which applies, in order, - the following sequence of transformations to the design weights: - (1) a bit mask for the boundary pixels, (2) convolution with a - conic filter, and (3) projection via a hyperbolic tangent (if - necessary). +def filter_projection( + weights: np.ndarray, sigmoid_threshold: float, sigmoid_bias: float +) -> np.ndarray: + """A differentiable function to filter and project the design weights. Args: - x: design weights as a 1d array of size Nx*Ny. - eta: erosion/dilation parameter for the projection. - beta: bias parameter for the projection. A value of 0 is no projection. + weights: design weights as a flattened (1D) array. + sigmoid_threshold: erosion/dilation parameter for the projection. + sigmoid_bias: bias parameter for the projection. 0 is no projection. Returns: - The mapped design weights as a 1d array. + The mapped design weights as a 1D array. """ - x = npa.where( - Si_mask.flatten(), + silicon_mask, silicon_dioxide_mask = border_masks() + + weights_masked = npa.where( + silicon_mask.flatten(), 1, npa.where( - SiO2_mask.flatten(), + silicon_dioxide_mask.flatten(), 0, - x, + weights, ), ) - filtered_field = mpa.conic_filter( - x, - filter_radius, - design_region_size.x, - design_region_size.y, - design_region_resolution, + weights_filtered = mpa.conic_filter( + weights_masked, + filter_radius_um, + DESIGN_REGION_UM.x, + DESIGN_REGION_UM.y, + DESIGN_REGION_RESOLUTION_UM, ) - if beta == 0: - return filtered_field.flatten() - + if sigmoid_bias == 0: + return weights_filtered.flatten() else: - projected_field = mpa.tanh_projection( - filtered_field, - beta, - eta, + weights_projected = mpa.tanh_projection( + weights_filtered, + sigmoid_bias, + sigmoid_threshold, ) + return weights_projected.flatten() - return projected_field.flatten() - -def f(x: np.ndarray, grad: np.ndarray) -> float: +def obj_func(epigraph_and_weights: np.ndarray, grad: np.ndarray) -> float: """Objective function for the epigraph formulation. Args: - x: 1d array of size 1+Nx*Ny containing epigraph variable (first element) - and design weights (remaining Nx*Ny elements). - grad: the gradient as a 1d array of size 1+Nx*Ny modified in place. + epigraph_and_weights: 1D array containing epigraph variable (first + element) and design weights (remaining elements). + grad: the gradient as a flattened (1D) array, modified in place. Returns: - The epigraph variable (a scalar). + The scalar epigraph variable. """ - t = x[0] # epigraph variable - v = x[1:] # design weights + epigraph = epigraph_and_weights[0] + if grad.size > 0: grad[0] = 1 grad[1:] = 0 - return t + + return epigraph -def c( +def epigraph_constraint( result: np.ndarray, - x: np.ndarray, + epigraph_and_weights: np.ndarray, gradient: np.ndarray, - eta: float, - beta: float, + sigmoid_threshold: float, + sigmoid_bias: float, use_epsavg: bool, -): +) -> None: """Constraint function for the epigraph formulation. Args: - result: the result of the function evaluation modified in place. - x: 1d array of size 1+Nx*Ny containing epigraph variable (first - element) and design weights (remaining Nx*Ny elements). - gradient: the Jacobian matrix with dimensions (1+Nx*Ny, - 2*num. wavelengths) modified in place. - eta: erosion/dilation parameter for projection. - beta: bias parameter for projection. + result: the result of the function evaluation, modified in place. + epigraph_and_weights: 1D array containing the epigraph variable (first + element) and design weights (remaining elements). + gradient: the Jacobian matrix with dimensions (1 + NX_DESIGN_GRID * + NY_DESIGN_GRID, 2 * num. wavelengths), modified in place. + sigmoid_threshold: erosion/dilation parameter for projection. + sigmoid_bias: bias parameter for projection. use_epsavg: whether to use subpixel smoothing. """ - t = x[0] # epigraph variable - v = x[1:] # design weights + epigraph = epigraph_and_weights[0] + weights = epigraph_and_weights[1:] - f0, dJ_du = opt([mapping(v, eta, 0 if use_epsavg else beta)]) + obj_val, grad = opt( + [ + filter_projection( + weights, sigmoid_threshold, 0 if use_epsavg else sigmoid_bias + ) + ] + ) - f0_reflection = f0[0] - f0_transmission = f0[1] - f0_merged = np.concatenate((f0_reflection, f0_transmission)) - f0_merged_str = "[" + ",".join(str(ff) for ff in f0_merged) + "]" + reflectance = obj_val[0] + transmittance = obj_val[1] + obj_val_merged = np.concatenate((reflectance, transmittance)) + obj_val_merged_str = str_from_list(obj_val_merged) - dJ_du_reflection = dJ_du[0] - dJ_du_transmission = dJ_du[1] - nfrq = len(frqs) - my_grad = np.zeros((Nx * Ny, 2 * nfrq)) - my_grad[:, :nfrq] = dJ_du_reflection - my_grad[:, nfrq:] = dJ_du_transmission + grad_reflectance = grad[0] + grad_transmittance = grad[1] + nfrq = len(frequencies) + grad = np.zeros((NX_DESIGN_GRID * NY_DESIGN_GRID, 2 * nfrq)) + grad[:, :nfrq] = grad_reflectance + grad[:, nfrq:] = grad_transmittance - # backpropagate the gradients through mapping function + # backpropagate the gradients through filter and projection function for k in range(2 * nfrq): - my_grad[:, k] = tensor_jacobian_product(mapping, 0)( - v, - eta, - beta, - my_grad[:, k], + grad[:, k] = tensor_jacobian_product(filter_projection, 0)( + weights, + sigmoid_threshold, + sigmoid_bias, + grad[:, k], ) if gradient.size > 0: - gradient[:, 0] = -1 # gradient w.r.t. epigraph variable ("t") - gradient[:, 1:] = my_grad.T # gradient w.r.t. each frequency objective + gradient[:, 0] = -1 # gradient w.r.t. epigraph variable + gradient[:, 1:] = grad.T # gradient w.r.t. each frequency objective - result[:] = np.real(f0_merged) - t + result[:] = np.real(obj_val_merged) - epigraph - objfunc_history.append(np.real(f0_merged)) - epivar_history.append(t) + objfunc_history.append(np.real(obj_val_merged)) + epivar_history.append(epigraph) print( - f"iteration:, {cur_iter[0]:3d}, eta: {eta}, beta: {beta:2d}, " - f"t: {t:.5f}, obj. func.: {f0_merged_str}" + f"iteration:, {cur_iter[0]:3d}, sigmoid_bias: {sigmoid_bias:2d}, " + f"epigraph: {epigraph:.5f}, obj. func.: {obj_val_merged_str}" ) cur_iter[0] = cur_iter[0] + 1 -def glc(result: np.ndarray, x: np.ndarray, gradient: np.ndarray, beta: float) -> float: +def line_width_and_spacing_constraint( + result: np.ndarray, + epigraph_and_weights: np.ndarray, + gradient: np.ndarray, + sigmoid_bias: float, +) -> float: """Constraint function for the minimum linewidth. Args: result: the result of the function evaluation modified in place. - x: 1d array of size 1+Nx*Ny containing epigraph variable (first - element) and design weights (remaining elements). - gradient: the Jacobian matrix with dimensions (1+Nx*Ny, - num. wavelengths) modified in place. - beta: bias parameter for projection. + epigraph_and_weights: 1D array containing epigraph variable (first + element) and design weights (remaining elements). + gradient: the Jacobian matrix, modified in place. + sigmoid_bias: bias parameter for projection. Returns: The value of the constraint function (a scalar). """ - t = x[0] # dummy parameter - v = x[1:] # design parameters + epigraph = epigraph_and_weights[0] + weights = epigraph_and_weights[1:] a1 = 1e-3 # hyper parameter (primary) b1 = 0 # hyper parameter (secondary) gradient[:, 0] = -a1 - filter_f = lambda a: mpa.conic_filter( - a.reshape(Nx, Ny), - filter_radius, - design_region_size.x, - design_region_size.y, - design_region_resolution, + filter_func = lambda a: mpa.conic_filter( + a.reshape(NX_DESIGN_GRID, NY_DESIGN_GRID), + filter_radius_um, + DESIGN_REGION_UM.x, + DESIGN_REGION_UM.y, + DESIGN_REGION_RESOLUTION_UM, + ) + + threshold_func = lambda a: mpa.tanh_projection( + a, sigmoid_bias, SIGMOID_THRESHOLD_INTRINSIC ) - threshold_f = lambda a: mpa.tanh_projection(a, beta, eta_i) # hyper parameter (constant factor and exponent) - c0 = 1e7 * (filter_radius * 1 / resolution) ** 4 + c0 = 1e7 * (filter_radius_um * 1 / RESOLUTION_UM) ** 4 - M1 = lambda a: mpa.constraint_solid(a, c0, eta_e, filter_f, threshold_f, 1) - M2 = lambda a: mpa.constraint_void(a, c0, eta_d, filter_f, threshold_f, 1) + M1 = lambda a: mpa.constraint_solid( + a, c0, SIGMOID_THRESHOLD_EROSION, filter_func, threshold_func, 1 + ) + M2 = lambda a: mpa.constraint_void( + a, c0, SIGMOID_THRESHOLD_DILATION, filter_func, threshold_func, 1 + ) - g1 = grad(M1)(v) - g2 = grad(M2)(v) + g1 = grad(M1)(weights) + g2 = grad(M2)(weights) - result[0] = M1(v) - a1 * t - b1 - result[1] = M2(v) - a1 * t - b1 + result[0] = M1(weights) - a1 * epigraph - b1 + result[1] = M2(weights) - a1 * epigraph - b1 gradient[0, 1:] = g1.flatten() gradient[1, 1:] = g2.flatten() - t1 = (M1(v) - b1) / a1 - t2 = (M2(v) - b1) / a1 + t1 = (M1(weights) - b1) / a1 + t2 = (M2(weights) - b1) / a1 - print(f"glc:, {result[0]}, {result[1]}, {t1}, {t2}") + print(f"line_width_and_spacing_constraint:, {result[0]}, {result[1]}, {t1}, {t2}") return max(t1, t2) def straight_waveguide() -> (np.ndarray, NamedTuple): - """Computes the DFT fields from the mode source in a straight waveguide - for use as normalization of the reflectance measurement during the - optimization. + """Computes the DFT fields from the mode source in a straight waveguide. + + The DFT fields are used as normalization of the reflectance measurement + during the optimization. Returns: A 2-tuple consisting of a 1d array of DFT fields and DFT fields object @@ -304,26 +316,26 @@ def straight_waveguide() -> (np.ndarray, NamedTuple): """ sources = [ mp.EigenModeSource( - src=mp.GaussianSource(fcen, fwidth=df), - size=mp.Vector3(0, sy, 0), + src=mp.GaussianSource(frequency_center, fwidth=frequency_width), + size=mp.Vector3(0, cell_um.y, 0), center=src_pt, eig_band=1, - eig_parity=eig_parity, + eig_parity=MODE_SYMMETRY, ) ] geometry = [ mp.Block( - size=mp.Vector3(mp.inf, w, mp.inf), + size=mp.Vector3(mp.inf, WAVEGUIDE_UM.y, mp.inf), center=mp.Vector3(), - material=Si, + material=SILICON, ) ] sim = mp.Simulation( - resolution=resolution, - default_material=SiO2, - cell_size=cell_size, + resolution=RESOLUTION_UM, + default_material=SILICON_DIOXIDE, + cell_size=cell_um, sources=sources, geometry=geometry, boundary_layers=pml_layers, @@ -331,8 +343,8 @@ def straight_waveguide() -> (np.ndarray, NamedTuple): ) refl_mon = sim.add_mode_monitor( - frqs, - mp.ModeRegion(center=refl_pt, size=mp.Vector3(0, sy, 0)), + frequencies, + mp.ModeRegion(center=refl_pt, size=mp.Vector3(0, cell_um.y, 0)), yee_grid=True, ) @@ -341,7 +353,7 @@ def straight_waveguide() -> (np.ndarray, NamedTuple): res = sim.get_eigenmode_coefficients( refl_mon, [1], - eig_parity=eig_parity, + eig_parity=MODE_SYMMETRY, ) coeffs = res.alpha @@ -356,12 +368,12 @@ def mode_converter_optimization( input_flux_data: NamedTuple, use_damping: bool, use_epsavg: bool, - beta: float, + sigmoid_bias: float, ) -> mpa.OptimizationProblem: """Sets up the adjoint optimization of the waveguide mode converter. Args: - input_flux: 1d array of DFT fields from normalization run. + input_flux: 1D array of DFT fields from the normalization run. input_flux_data: DFT fields object returned by `meep.get_flux_data`. use_damping: whether to use the damping feature of `MaterialGrid`. use_epsavg: whether to use subpixel smoothing in `MaterialGrid`. @@ -370,21 +382,18 @@ def mode_converter_optimization( A `meep.adjoint.OptimizationProblem` class object. """ matgrid = mp.MaterialGrid( - mp.Vector3(Nx, Ny, 0), - SiO2, - Si, - weights=np.ones((Nx, Ny)), - beta=beta if use_epsavg else 0, + mp.Vector3(NX_DESIGN_GRID, NY_DESIGN_GRID, 0), + SILICON_DIOXIDE, + SILICON, + weights=np.ones((NX_DESIGN_GRID, NY_DESIGN_GRID)), + beta=sigmoid_bias if use_epsavg else 0, do_averaging=True if use_epsavg else False, - damping=0.02 * 2 * np.pi * fcen if use_damping else 0, + damping=0.02 * 2 * np.pi * frequency_center if use_damping else 0, ) matgrid_region = mpa.DesignRegion( matgrid, - volume=mp.Volume( - center=mp.Vector3(), - size=mp.Vector3(design_region_size.x, design_region_size.y, mp.inf), - ), + volume=mp.Volume(center=mp.Vector3(), size=DESIGN_REGION_UM), ) matgrid_geometry = [ @@ -398,8 +407,8 @@ def mode_converter_optimization( geometry = [ mp.Block( center=mp.Vector3(), - size=mp.Vector3(mp.inf, w, mp.inf), - material=Si, + size=mp.Vector3(mp.inf, WAVEGUIDE_UM.y, mp.inf), + material=SILICON, ) ] @@ -407,18 +416,18 @@ def mode_converter_optimization( sources = [ mp.EigenModeSource( - src=mp.GaussianSource(fcen, fwidth=df), - size=mp.Vector3(0, sy, 0), + src=mp.GaussianSource(frequency_center, fwidth=frequency_width), + size=mp.Vector3(0, cell_um.y, 0), center=src_pt, eig_band=1, - eig_parity=eig_parity, + eig_parity=MODE_SYMMETRY, ), ] sim = mp.Simulation( - resolution=resolution, - default_material=SiO2, - cell_size=cell_size, + resolution=RESOLUTION_UM, + default_material=SILICON_DIOXIDE, + cell_size=cell_um, sources=sources, geometry=geometry, boundary_layers=pml_layers, @@ -430,21 +439,21 @@ def mode_converter_optimization( sim, mp.Volume( center=refl_pt, - size=mp.Vector3(0, sy, 0), + size=mp.Vector3(0, cell_um.y, 0), ), 1, forward=False, - eig_parity=eig_parity, + eig_parity=MODE_SYMMETRY, subtracted_dft_fields=input_flux_data, ), mpa.EigenmodeCoefficient( sim, mp.Volume( center=tran_pt, - size=mp.Vector3(0, sy, 0), + size=mp.Vector3(0, cell_um.y, 0), ), 2, - eig_parity=eig_parity, + eig_parity=MODE_SYMMETRY, ), ] @@ -459,7 +468,7 @@ def J2(refl_mon, tran_mon): objective_functions=[J1, J2], objective_arguments=obj_list, design_regions=[matgrid_region], - frequencies=frqs, + frequencies=frequencies, ) return opt @@ -468,114 +477,132 @@ def J2(refl_mon, tran_mon): if __name__ == "__main__": input_flux, input_flux_data = straight_waveguide() - algorithm = nlopt.LD_MMA + silicon_mask, silicon_dioxide_mask = border_masks() - # number of design parameters - n = Nx * Ny + num_weights = NX_DESIGN_GRID * NY_DESIGN_GRID - # initial design parameters - x = np.ones((n,)) * 0.5 - x[Si_mask.flatten()] = 1.0 # set the edges of waveguides to silicon - x[SiO2_mask.flatten()] = 0.0 # set the other edges to SiO2 + # Initial design weights. + epigraph_and_weights = np.ones((num_weights,)) * 0.5 + epigraph_and_weights[silicon_mask.flatten()] = 1.0 + epigraph_and_weights[silicon_dioxide_mask.flatten()] = 0.0 - # lower and upper bounds for design weights - lb = np.zeros((n,)) - lb[Si_mask.flatten()] = 1.0 - ub = np.ones((n,)) - ub[SiO2_mask.flatten()] = 0.0 + # Lower and upper bounds for design weights. + weights_lower_bound = np.zeros((num_weights,)) + weights_lower_bound[silicon_mask.flatten()] = 1.0 + weights_upper_bound = np.ones((num_weights,)) + weights_upper_bound[silicon_dioxide_mask.flatten()] = 0.0 - # insert epigraph variable initial value (arbitrary) and bounds into the - # design array. the actual value is determined by the objective and + # Insert epigraph variable with arbitrary initial value and bounds into the + # design array. The actual value is determined by the objective and # constraint functions below. - x = np.insert(x, 0, 1.2) - lb = np.insert(lb, 0, -np.inf) - ub = np.insert(ub, 0, +np.inf) + epigraph_and_weights = np.insert(epigraph_and_weights, 0, 1.2) + weights_lower_bound = np.insert(weights_lower_bound, 0, -np.inf) + weights_upper_bound = np.insert(weights_upper_bound, 0, +np.inf) objfunc_history = [] epivar_history = [] cur_iter = [0] - beta_thresh = 64 # threshold beta above which to use subpixel smoothing - betas = [8, 16, 32, 64, 128, 256] + sigmoid_bias_threshold = 64 # threshold beta above which to use subpixel smoothing + sigmoid_biases = [8, 16, 32, 64, 128, 256] max_evals = [80, 80, 100, 120, 120, 100] - tol_epi = np.array([1e-4] * 2 * len(frqs)) # R, 1-T - tol_lw = np.array([1e-8] * 2) # line width, line spacing - - for beta, max_eval in zip(betas, max_evals): - solver = nlopt.opt(algorithm, n + 1) - solver.set_lower_bounds(lb) - solver.set_upper_bounds(ub) - solver.set_min_objective(f) + tolerance_epigraph = np.array([1e-4] * 2 * len(frequencies)) # R, 1-T + tolerance_width_and_spacing = np.array([1e-8] * 2) # line width, line spacing + + for sigmoid_bias, max_eval in zip(sigmoid_biases, max_evals): + solver = nlopt.opt(nlopt.LD_MMA, num_weights + 1) + solver.set_lower_bounds(weights_lower_bound) + solver.set_upper_bounds(weights_upper_bound) + solver.set_min_objective(obj_func) solver.set_maxeval(max_eval) solver.set_param("dual_ftol_rel", 1e-7) solver.add_inequality_mconstraint( - lambda rr, xx, gg: c( - rr, - xx, - gg, - eta_i, - beta, - False if beta < beta_thresh else True, + lambda result_, epigraph_and_weights_, grad_: epigraph_constraint( + result_, + epigraph_and_weights_, + grad_, + SIGMOID_THRESHOLD_INTRINSIC, + sigmoid_bias, + False if sigmoid_bias < sigmoid_bias_threshold else True, ), - tol_epi, + tolerance_epigraph, ) solver.set_param("verbosity", 1) + if sigmoid_bias < sigmoid_bias_threshold: + use_epsavg = False + else: + use_epsavg = True + opt = mode_converter_optimization( input_flux, input_flux_data, True, # use_damping - False if beta < beta_thresh else True, # use_epsavg - beta, + use_epsavg, + sigmoid_bias, ) - # apply the minimum linewidth constraint + # Apply the minimum linewidth constraint # only in the final epoch to an initial # binary design from the previous epoch. - if beta == betas[-1]: - res = np.zeros(2) - grd = np.zeros((2, n + 1)) - t = glc(res, x, grd, beta) + if sigmoid_bias == sigmoid_biases[-1]: + line_width_and_spacing = np.zeros(2) + grad_line_width_and_spacing = np.zeros((2, num_weights + 1)) + linewidth_constraint_val = line_width_and_spacing_constraint( + line_width_and_spacing, + epigraph_and_weights, + grad_line_width_and_spacing, + sigmoid_bias, + ) solver.add_inequality_mconstraint( - lambda rr, xx, gg: glc( - rr, - xx, - gg, - beta, + lambda result_, epigraph_and_weights_, grad_: line_width_and_spacing_constraint( + result_, + epigraph_and_weights_, + grad_, + sigmoid_bias, ), - tol_lw, + tolerance_width_and_spacing, ) - # execute a single forward run before the start of each + # Execute a single forward run before the start of each # epoch and manually set the initial epigraph variable to # slightly larger than the largest value of the objective # function over the six wavelengths and the lengthscale # constraint (final epoch only). - t0 = opt( + epigraph_initial = opt( [ - mapping( - x[1:], - eta_i, - beta if beta < beta_thresh else 0, + filter_projection( + epigraph_and_weights[1:], + SIGMOID_THRESHOLD_INTRINSIC, + sigmoid_bias if sigmoid_bias < sigmoid_bias_threshold else 0, ), ], need_gradient=False, ) - t0 = np.concatenate((t0[0][0], t0[0][1])) - t0_str = "[" + ",".join(str(tt) for tt in t0) + "]" - x[0] = np.amax(t0) - x[0] = 1.05 * (max(x[0], t) if beta == betas[-1] else x[0]) - print(f"data:, {beta}, {t0_str}, {x[0]}") + epigraph_initial = np.concatenate( + (epigraph_initial[0][0], epigraph_initial[0][1]) + ) + epigraph_initial_str = str_from_list(epigraph_initial) + + epigraph_and_weights[0] = np.amax(epigraph_initial) + if sigmoid_bias == sigmoid_biases[-1]: + epigraph_and_weights[0] = 1.05 * max( + epigraph_and_weights[0], linewidth_constraint_val + ) + print( + f"epigraph-calibration:, {sigmoid_bias}, {epigraph_initial_str}, " + f"{epigraph_and_weights[0]}" + ) - x[:] = solver.optimize(x) + epigraph_and_weights[:] = solver.optimize(epigraph_and_weights) - optimal_design_weights = mapping( - x[1:], - eta_i, - beta, - ).reshape(Nx, Ny) + optimal_design_weights = filter_projection( + epigraph_and_weights[1:], + SIGMOID_THRESHOLD_INTRINSIC, + sigmoid_bias, + ).reshape(NX_DESIGN_GRID, NY_DESIGN_GRID) - # save the unmapped weights and a bitmap image + # Save the unmapped weights and a bitmap image # of the design weights at the end of each epoch. fig, ax = plt.subplots() ax.imshow( @@ -586,33 +613,38 @@ def J2(refl_mon, tran_mon): ax.set_axis_off() if mp.am_master(): fig.savefig( - f"optimal_design_beta{beta}.png", + f"optimal_design_beta{sigmoid_bias}.png", dpi=150, bbox_inches="tight", ) - # save the final (unmapped) design as a 2d array in CSV format + # Save the final (unmapped) design as a 2D array in CSV format np.savetxt( - f"unmapped_design_weights_beta{beta}.csv", - x[1:].reshape(Nx, Ny), + f"unmapped_design_weights_beta{sigmoid_bias}.csv", + epigraph_and_weights[1:].reshape(NX_DESIGN_GRID, NY_DESIGN_GRID), fmt="%4.2f", delimiter=",", ) - # save all the important optimization parameters and output - # as separate arrays in a single file for post processing. - with open("optimal_design.npz", "wb") as fl: + # Save important optimization parameters and output as + # separate fields for post processing. + with open("optimal_design.npz", "wb") as datafile: np.savez( - fl, - Nx=Nx, - Ny=Ny, - design_region_size=(dx, dy), - design_region_resolution=design_region_resolution, - betas=betas, + datafile, + RESOLUTION_UM=RESOLUTION_UM, + DESIGN_WAVELENGTHS_UM=DESIGN_WAVELENGTHS_UM, + WAVEGUIDE_UM=WAVEGUIDE_UM, + PADDING_UM=PADDING_UM, + PML_UM=PML_UM, + DESIGN_REGION_UM=DESIGN_REGION_UM, + DESIGN_REGION_RESOLUTION_UM=DESIGN_REGION_RESOLUTION_UM, + NX_DESIGN_GRID=NX_DESIGN_GRID, + NY_DESIGN_GRID=NY_DESIGN_GRID, + MIN_LENGTH_UM=MIN_LENGTH_UM, + sigmoid_biases=sigmoid_biases, max_eval=max_eval, objfunc_history=objfunc_history, epivar_history=epivar_history, - t=x[0], - unmapped_design_weights=x[1:], - minimum_length=minimum_length, + epigraph_variable=epigraph_and_weights[0], + unmapped_design_weights=epigraph_and_weights[1:], optimal_design_weights=optimal_design_weights, )