diff --git a/metapeg.ipynb b/metapeg.ipynb new file mode 100644 index 0000000..8e88383 --- /dev/null +++ b/metapeg.ipynb @@ -0,0 +1,669 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "8658f8cd", + "metadata": {}, + "source": [ + "This code implements an inverse designed metasurface that is used to recreate a given intenstity\n", + "pattern. It is based on meep and tidy3d tutorials." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4c8450d1", + "metadata": {}, + "outputs": [], + "source": [ + "import meep as mp\n", + "import meep.adjoint as mpa\n", + "import numpy as np\n", + "from autograd import numpy as npa\n", + "from autograd import tensor_jacobian_product, grad\n", + "import nlopt\n", + "from matplotlib import pyplot as plt\n", + "from matplotlib.patches import Circle\n", + "import scipy as sp\n", + "import os\n", + "\n", + "mp.verbosity(0)\n", + "\n", + "#setup a results folder\n", + "if not os.path.exists(\"results\"):\n", + " os.makedirs(\"results\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c728065b", + "metadata": {}, + "outputs": [], + "source": [ + "######PHYSICS SETUP######\n", + "subs_index = 3.4\n", + "subs_perm = subs_index**2\n", + "\n", + "Si = mp.Medium(index = subs_index)\n", + "Air = mp.Medium(index = 1.0)\n", + "\n", + "wavelength = 1.55\n", + "nf = 1 #number of frequencies to test\n", + "frequencies = np.array([1 / 1.55])\n", + "\n", + "#distance between PML and source/monitor\n", + "buffer = 2.5\n", + "\n", + "######PML######\n", + "pml_size = 1.0\n", + "pml_layers = [mp.PML(pml_size)]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ca9ab89e", + "metadata": {}, + "outputs": [], + "source": [ + "######FILTER######\n", + "minimum_length = 0.09 # minimum length scale (microns)\n", + "eta_i = (\n", + " 0.5 # blueprint (or intermediate) design field thresholding point (between 0 and 1)\n", + ")\n", + "eta_e = 0.55 # erosion design field thresholding point (between 0 and 1)\n", + "eta_d = 1 - eta_e # dilation design field thresholding point (between 0 and 1)\n", + "filter_radius = mpa.get_conic_radius_from_eta_e(minimum_length, eta_e)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5fc04b92", + "metadata": {}, + "outputs": [], + "source": [ + "######DESIGN REGION######\n", + "design_region_width = 3.0\n", + "design_region_length = 3.0\n", + "design_region_height = (wavelength/np.sqrt(subs_index))/2 #half-wavelength in the substrate\n", + "\n", + "resolution = 20 # pixels/um\n", + "design_region_resolution = int(resolution)\n", + "\n", + "#number of pixels in the design region\n", + "Nx = int(design_region_width * design_region_resolution) + 1\n", + "Ny = int(design_region_length * design_region_resolution) + 1\n", + "\n", + "design_variables = mp.MaterialGrid(mp.Vector3(Nx, Ny), Air, Si, grid_type=\"U_MEAN\")\n", + "\n", + "design_region = mpa.DesignRegion(\n", + " design_variables,\n", + " volume=mp.Volume(\n", + " center = mp.Vector3(0,0,design_region_height / 2),\n", + " size = mp.Vector3(\n", + " design_region_width, \n", + " design_region_length, \n", + " design_region_height),\n", + " ),\n", + ")\n", + "\n", + "######CELL GEOMETRY######\n", + "Sx = 2 * pml_size + design_region_width\n", + "Sy = 2 * pml_size + design_region_length\n", + "Sz = 2 * pml_size + design_region_height + 4*buffer\n", + "cell_size = mp.Vector3(Sx, Sy, Sz)\n", + "\n", + "substrate_thickness = Sz / 2\n", + "substrate_size = mp.Vector3(design_region_width + 2 * pml_size, \n", + " design_region_length + 2 * pml_size, \n", + " substrate_thickness)\n", + "\n", + "substrate_center = mp.Vector3(0, 0, - substrate_thickness / 2)\n", + "\n", + "substrate_block = mp.Block(\n", + " center = substrate_center,\n", + " size = substrate_size,\n", + " material = Si,\n", + ")\n", + "\n", + "design_region_block = mp.Block(\n", + " center = design_region.center, \n", + " size = design_region.size, \n", + " material = design_variables,\n", + ")\n", + "\n", + "geometry = [#combine geometries\n", + " substrate_block,\n", + " design_region_block,\n", + "]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3163836f", + "metadata": {}, + "outputs": [], + "source": [ + "######SOURCE######\n", + "fcen = 1 / wavelength\n", + "width = 0.2\n", + "fwidth = width * fcen\n", + "source_center = [0, 0, -3]\n", + "source_size = mp.Vector3(design_region_width, design_region_length, 0)\n", + "src = mp.GaussianSource(frequency=fcen, fwidth=fwidth)\n", + "source = [mp.Source(src, component=mp.Ex, size=source_size, center=source_center)]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d2f566e1", + "metadata": {}, + "outputs": [], + "source": [ + "######OBSERVATION PLANES######\n", + "#planes used for references in 2D plots\n", + "output_plane = mp.Volume(\n", + " center=mp.Vector3(0,0,0), \n", + " size=mp.Vector3(Sx,0,Sz))\n", + "\n", + "output_plane_xy = mp.Volume(\n", + " center=mp.Vector3(0,0,design_region_height / 2),\n", + " size=mp.Vector3(Sx,Sy,0)\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "41b6bdaa", + "metadata": {}, + "outputs": [], + "source": [ + "######SIMULATION######\n", + "kpoint = mp.Vector3()\n", + "sim = mp.Simulation(\n", + " cell_size = cell_size,\n", + " boundary_layers = pml_layers,\n", + " geometry = geometry,\n", + " sources = source,\n", + " default_material = Air,\n", + " resolution = resolution,\n", + ")\n", + "\n", + "sim.plot2D(output_plane=output_plane,show_sources=True, show_boundaries=True, show_geometry=True)\n", + "plt.savefig(\"results/simulation_geometry.pdf\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3c05d8e1", + "metadata": {}, + "outputs": [], + "source": [ + "import xarray as xr\n", + "from PIL import Image\n", + "\n", + "logo_name = \"misc/peg.png\"\n", + "\n", + "def get_logo() -> np.ndarray:\n", + " logo = Image.open(logo_name).convert(\"L\") #L mode stands for grayscale\n", + " logo = np.array(logo).astype(float)\n", + " #logo is normalized between 0 and 1\n", + " logo -= np.min(logo)\n", + " logo /= np.max(logo)\n", + " #in general, logo has any shape. At a later stage it will\n", + " #be rescaled to fit the design region\n", + " print(\"Logo shape: \", logo.shape)\n", + " return logo\n", + "\n", + "def intensity_desired_fn_logo(xs: list, ys: list, rescale: float = 0.5) -> np.ndarray:\n", + " '''return the value of the logo as a function of the x,y coordinates, \n", + " with a rescale factor\n", + " Parameters:\n", + " xs: list of x coordinates\n", + " ys: list of y coordinates\n", + " Returns:\n", + " logo_interp: 2D array of the same shape as the design region (nx,ny)\n", + " '''\n", + " logo_values = get_logo()\n", + " logo_values = 1 - np.rot90(np.rot90(np.rot90(logo_values))) #rotate the logo by 90 degrees to match the orientation of the simulation\n", + " \n", + " nx,ny = logo_values.shape\n", + " xs_logo = np.linspace(rescale*min(xs), rescale*max(xs), nx)\n", + " ys_logo = np.linspace(rescale*min(ys), rescale*max(ys), ny)\n", + " logo_dataarray = xr.DataArray(logo_values, coords=dict(x=xs_logo, y=ys_logo))\n", + "\n", + " logo_interp = logo_dataarray.interp(x=xs, y=ys)\n", + "\n", + " return np.nan_to_num(logo_interp.values, nan=np.min(logo_interp))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "456d598e", + "metadata": {}, + "outputs": [], + "source": [ + "xs = ys = np.linspace(-design_region_width / 2, design_region_width / 2, Nx)\n", + "intensity_desired = intensity_desired_fn_logo(xs,ys)\n", + "print(intensity_desired.shape)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5c6ed0bc", + "metadata": {}, + "outputs": [], + "source": [ + "plt.pcolormesh(xs, ys, intensity_desired.T, cmap=\"magma\")\n", + "plt.gca().set_aspect(\"equal\")\n", + "plt.xlabel(\"x\")\n", + "plt.ylabel(\"y\")\n", + "plt.title(\"desired intensity pattern\")\n", + "plt.colorbar()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8d1d7dc0", + "metadata": {}, + "outputs": [], + "source": [ + "######MONITOR######\n", + "far_field_obs_point_z_coord = 10\n", + "\n", + "plane_points_array = np.array([[x_p,y_p] for x_p in xs for y_p in ys])\n", + "far_z_plane = [mp.Vector3(x_p,y_p,far_field_obs_point_z_coord) for x_p in xs for y_p in ys] \n", + "plane_points_matrix = plane_points_array.reshape(Nx,Ny,2)\n", + "#make a plot of the sampling coordinates\n", + "plt.scatter(plane_points_array[:,0], plane_points_array[:,1], s=1)\n", + "plt.xlabel(\"x\")\n", + "plt.ylabel(\"y\")\n", + "plt.title(\"sampling coordinates for near-to-far transformation\")\n", + "\n", + "#define the near-field observation plane (that spans the whole design region)\n", + "NearRegions = [mp.Near2FarRegion(\n", + " center=mp.Vector3(0,0, design_region_height + wavelength),\n", + " size=mp.Vector3(design_region_width,design_region_length,0),\n", + " weight=+1,\n", + ")]\n", + "\n", + "#mpa.Near2FarFields return a numpy array with shape (num_of_points, nfreq, 6) where the \n", + "#third axis are the field components Ex (0),Ey (1),Ez (2),Hx (3),Hy (4),Hz (5)\n", + "FarFields = mpa.Near2FarFields(sim, NearRegions, far_z_plane)\n", + "ob_list = [FarFields]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c725fa9c", + "metadata": {}, + "outputs": [], + "source": [ + "######INTENSITIES######\n", + "def get_intensities(FF:np.ndarray) -> tuple[np.ndarray, np.ndarray]:\n", + " '''Utility function to compute target and measured intensities given far field data.\n", + " Input is a numpy array with shape (num_of_points, nfreq, 6) where the \n", + " third axis are the field components Ex, Ey, Ez, Hx, Hy, Hz. num_of_points represents the number\n", + " of far-field points that were passed in the list far_z_plane when defining the Near2FarRegions.\n", + " \n", + " Parameters:\n", + " FF: np.ndarray with shape (num_of_points, nfreq, 6)\n", + " '''\n", + " #first, retrieve the target intensities from the target\n", + " target_intensities = intensity_desired_fn_logo(xs,ys)\n", + " #second, compute the measured intensities from the far-field\n", + " #TODO need to check if the measured matrix needs to be flipped up/down or\n", + " #other transformations due to how it is constructed.\n", + " measured_intensities = npa.reshape(npa.abs(FF[:,:,0])**2, (Nx,Ny))\n", + " \n", + " return target_intensities, measured_intensities" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "362842c5", + "metadata": {}, + "outputs": [], + "source": [ + "######NORMALIZATION######\n", + "#compute the initial far-field intensity without \n", + "#the presence of the metasurface to use as normalization\n", + "norm_sim = mp.Simulation(\n", + " cell_size = cell_size,\n", + " boundary_layers = pml_layers,\n", + " geometry = geometry,\n", + " sources = source,\n", + " default_material = Air,\n", + " resolution = resolution,\n", + ")\n", + "\n", + "norm_near2far = norm_sim.add_near2far(fcen,0,1,*NearRegions)\n", + "\n", + "norm_sim.plot2D(output_plane=output_plane,show_sources=True, show_boundaries=True, show_geometry=True)\n", + "plt.savefig(\"results/normalization_simulation_geometry.pdf\")\n", + "\n", + "norm_sim.run(until=200)\n", + "ref_fields = np.array([norm_sim.get_farfield(norm_near2far, point) for point in far_z_plane])\n", + "#TODO this todo is linked to the one above. Need to check if the far-field need to be flipped or something\n", + "#at this time, the two definition are consistent with each other.\n", + "intensity_norm_mean = npa.mean(npa.abs(ref_fields[:,:,0])**2)\n", + "print(\"intensity_norm_mean shape (should be 1)\",intensity_norm_mean.shape)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "710abef1", + "metadata": {}, + "outputs": [], + "source": [ + "#this the objective to function\n", + "#the objective is to minimize the distance between the measured intenstity\n", + "# and the target.\n", + "intensity_range = int_min, int_max = (0.0, 1.0)\n", + "\n", + "def J1(FF:np.ndarray):\n", + " '''Computes the objective function for the given input field\n", + " Parameters: \n", + " FF: np.ndarray with shape (num_of_points, nfreq, 6)\n", + " Returns:\n", + " the value of the objective function, a float\n", + " '''\n", + " print(\"FF shape\", FF.shape)\n", + " target_intensity, measured_intensity = get_intensities(FF)\n", + "\n", + " measured_intensity = measured_intensity * npa.mean(target_intensity) / intensity_norm_mean\n", + "\n", + " int_range_magnitude = npa.abs(int_max - int_min)\n", + " measured_intensity = npa.clip(measured_intensity, int_min, int_max)\n", + " target_intensity = int_range_magnitude * intensity_desired + int_min\n", + "\n", + " difference = measured_intensity - target_intensity\n", + " difference_denominator = npa.full(target_intensity.shape, int_range_magnitude)\n", + "\n", + " #obj_to_return = npa.mean(npa.abs(FF[0,:,0]) ** 2)\n", + "\n", + " return npa.linalg.norm(difference) / npa.linalg.norm(difference_denominator)\n", + "#returns the first point of all frequencies of the Ex component" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "aa694ac2", + "metadata": {}, + "outputs": [], + "source": [ + "######OPTIMIZATION PROBLEM######\n", + "'''\n", + "An OptimizationProblem class object knows how to do one basic thing: \n", + "Given an input vector of design variables, compute the objective function value \n", + "(forward calculation) and optionally its gradient (adjoint calculation)\n", + "The gradient of the objective function is computed against its arguments, which are\n", + "the FarFields.\n", + "'''\n", + "opt = mpa.OptimizationProblem(\n", + " simulation=sim, #the optimizer uses this simulation object\n", + " objective_functions=[J1], #the optimizer call this objective function\n", + " objective_arguments=ob_list, #the optimizer passes this list of arguments to the objective function\n", + " design_regions=[design_region], #the optimizer updates the design variables in this design region\n", + " frequencies=frequencies, #the optimizer runs the simulation at these frequencies\n", + " maximum_run_time=100, #maximum optimizer run time\n", + ")\n", + "######CHECK PLOT######\n", + "figure = opt.plot2D(True, output_plane=output_plane)\n", + "plt.savefig(\"results/simulation_setup.pdf\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c182593d", + "metadata": {}, + "outputs": [], + "source": [ + "def mapping(x,eta,beta):\n", + " '''\n", + " Arguments:\n", + " x: design variables (grey scale)\n", + " eta: thresholding point for projection (between 0 and 1)\n", + " beta: binarization strength for projection (positive)\n", + " Returns:\n", + " projected_field: the filtered and projected design variables\n", + " '''\n", + " #mapping is necessary to reach a physical design.\n", + " #we start from the design variables (grey scale), which are values in the range (0,1).\n", + " #typically, we set an initial value of 0.5.\n", + " #filtering makes a local point behave more like its neighbourhood, avoiding checker patterns\n", + " filtered_field = mpa.conic_filter(\n", + " x,\n", + " filter_radius,\n", + " design_region_width,\n", + " design_region_length,\n", + " [design_region_resolution],\n", + " )\n", + "\n", + " #projection\n", + " #this is an hyperbolic tangent projection\n", + " #after filtering, the design is still blurry, so we project it\n", + " #using an tanh projection method to push the values towards the upper/lower bound\n", + " #If a filtered value is >0.5, then we push it to 1, otherwise, if <0.5, we push it to 0.\n", + " #parameter beta represents the degree of binarization,\n", + " #parameter eta is threshold point\n", + " projected_field = mpa.tanh_projection(filtered_field,beta,eta)\n", + " \n", + " #the following lines introduce symmetry. Need to understand if necessary \n", + " #projected_field = (\n", + " # npa.flipud(projected_field) + projected_field\n", + " #)/2\n", + "\n", + " #interpolate to actual materials\n", + " return projected_field.flatten()\n", + "\n", + "evaluation_history = []\n", + "cur_iter = [0]\n", + "\n", + "def f(v, gradient, cur_beta):\n", + " print(\"current iteration: {}\".format(cur_iter[0]+1))\n", + "\n", + " #value and grad\n", + " f0, dJ_du = opt([mapping(v,eta_i,cur_beta)])\n", + " if gradient.size > 0:\n", + " # NOTE: gradient is assigned in-place! No need to return.\n", + " gradient[...] = tensor_jacobian_product(mapping, 0)(\n", + " v, eta_i, cur_beta, np.sum(dJ_du, axis=1)\n", + " )\n", + " \n", + " evaluation_history.append(np.real(f0))\n", + " \n", + " fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(18, 5))\n", + " \n", + " opt.plot2D(\n", + " False, \n", + " ax=ax1,\n", + " output_plane=output_plane,\n", + " show_sources = False,\n", + " show_monitors = False,\n", + " show_boundary_layers = False,\n", + " )\n", + " circ1 = Circle((2,2), minimum_length/2)\n", + " ax1.add_patch(circ1)\n", + " \n", + " opt.plot2D(\n", + " False, \n", + " ax=ax2,\n", + " output_plane=output_plane_xy,\n", + " show_sources = False,\n", + " show_monitors = False,\n", + " show_boundary_layers = False,\n", + " )\n", + " circ2 = Circle((2,2), minimum_length/2)\n", + " ax2.add_patch(circ2)\n", + " \n", + " ax3.plot(evaluation_history)\n", + " ax3.set_xlabel(\"iteration\")\n", + " ax3.set_ylabel(\"objective function value\")\n", + " ax3.set_title(\"objective function history\")\n", + "\n", + " plt.savefig(\"results/iter_{}.pdf\".format(cur_iter[0]+1))\n", + "\n", + " cur_iter[0] += 1\n", + "\n", + " return np.real(f0)\n", + "\n", + "\n", + "algorithm = nlopt.LD_MMA\n", + "n = Nx * Ny #number of parameters to optimize\n", + "\n", + "x = np.ones((n,)) * 0.5 #initial guess\n", + "#x = np.insert(x, 0, -1)\n", + "#lower and upper bounds for design variables\n", + "lb = np.zeros((Nx * Ny,))\n", + "ub = np.ones((Nx * Ny,))\n", + "\n", + "cur_beta = 4 #starting beta\n", + "beta_scale = 2 #beta scaling (update) factor\n", + "num_betas = 6 #number of betas to test\n", + "update_factor = 12 #number of iterations per beta\n", + "#total number of simulations = num_betas * update_factor * 2 (forward + adjoint)\n", + "ftol = 1e-5\n", + "\n", + "optimization_start = True\n", + "if optimization_start:\n", + " for iters in range(num_betas):\n", + " solver = nlopt.opt(algorithm, n)\n", + " solver.set_lower_bounds(lb)\n", + " solver.set_upper_bounds(ub)\n", + " solver.set_min_objective(lambda v, g: f(v, g, cur_beta))\n", + " solver.set_maxeval(update_factor) #number of iterations for each beta\n", + " solver.set_ftol_rel(ftol)\n", + " #observe that x is updated in place, so at\n", + " # the end it will contain the optimized parameters\n", + " x[...] = solver.optimize(x)\n", + " cur_beta = cur_beta * beta_scale #update beta\n", + " print(\"current beta: \", cur_beta)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5225b3d4", + "metadata": {}, + "outputs": [], + "source": [ + "######RESULTS######\n", + "#plot the figure of merit over iterations\n", + "plt.figure()\n", + "plt.plot(evaluation_history)\n", + "plt.xlabel(\"Iteration\")\n", + "plt.ylabel(\"Figure of Merit\")\n", + "plt.grid()\n", + "plt.savefig(\"results/figure_of_merit.pdf\")\n", + "\n", + "#plot and save the final design\n", + "np.savetxt(\"results/final_design.txt\", mapping(x, eta_i, cur_beta))\n", + "opt.update_design([mapping(x, eta_i, cur_beta)])\n", + "\n", + "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))\n", + " \n", + "opt.plot2D(\n", + " False, \n", + " ax=ax1,\n", + " output_plane=output_plane,\n", + " show_sources = False,\n", + " show_monitors = False,\n", + " show_boundary_layers = False,\n", + ")\n", + "\n", + "circ1 = Circle((2,2), minimum_length/2)\n", + "ax1.add_patch(circ1)\n", + "ax1.axis(\"off\")\n", + "\n", + "opt.plot2D(\n", + " False, \n", + " ax=ax2,\n", + " output_plane=output_plane_xy,\n", + " show_sources = False,\n", + " show_monitors = False,\n", + " show_boundary_layers = False,\n", + ")\n", + "circ2 = Circle((2,2), minimum_length/2)\n", + "ax2.add_patch(circ2)\n", + "ax2.axis(\"off\")\n", + " \n", + "plt.savefig(\"results/final_design_xy.pdf\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5c0ec83f", + "metadata": {}, + "outputs": [], + "source": [ + "######FULL SIMULATION######\n", + "#after completing the optimization we define the full domain and \n", + "#verify that the intensity at the desired location matches the target.\n", + "#Recall that the far field plane was defined at a distance d=10\n", + "opt.update_design([mapping(x, eta_i, cur_beta)])\n", + "\n", + "full_Sz = 26\n", + "\n", + "opt.sim = mp.Simulation(\n", + " cell_size = mp.Vector3(Sx,Sy,full_Sz),\n", + " boundary_layers = pml_layers,\n", + " k_point = mp.Vector3(),\n", + " geometry = geometry,\n", + " sources = source,\n", + " default_material = Air,\n", + " resolution = resolution,\n", + ")\n", + "\n", + "src = mp.ContinuousSource(frequency=fcen, fwidth=fwidth)\n", + "source = [mp.Source(src, component=mp.Ex, size=source_size, center=source_center)]\n", + "opt.sim.change_sources(source)\n", + "\n", + "full_output_plane = mp.Volume(\n", + " center=mp.Vector3(0,0,0), \n", + " size=mp.Vector3(Sx,0,full_Sz))\n", + "\n", + "opt.sim.plot2D(output_plane=full_output_plane,fields=mp.Ex)\n", + "plt.savefig(\"results/full_simulation_geometry.pdf\")\n", + "opt.sim.run(until=100)\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "mp", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.13.12" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}