commit 005b69cec47253f6c0eb3f1e0bbd422645c4f511 Author: edoardo Date: Tue May 12 23:13:08 2026 +0200 first commit diff --git a/mode_converter.ipynb b/mode_converter.ipynb new file mode 100644 index 0000000..f20f76c --- /dev/null +++ b/mode_converter.ipynb @@ -0,0 +1,459 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "f65ac12a", + "metadata": {}, + "source": [ + "This is a first attempt at creating an inverse designed fiber mode converter in meep. At first, we try to implement and combine tutorials from both tidy3d and meep." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7b374471", + "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" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3955252a", + "metadata": {}, + "outputs": [], + "source": [ + "######PHYSICS SETUP######\n", + "\n", + "mp.verbosity(1)\n", + "\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 = 3 #number of frequencies to test\n", + "frequencies = np.array([1/1.5, 1 / 1.55, 1 / 1.6])\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": "b39d1f98", + "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": "b1317dd3", + "metadata": {}, + "outputs": [], + "source": [ + "######DESIGN REGION######\n", + "design_region_width = 4.0\n", + "design_region_length = 4.0\n", + "design_region_height = 1.0\n", + "\n", + "resolution = 15 # 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(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 + 2*buffer\n", + "print(\"Cell size: \", Sx, Sy, Sz)\n", + "cell_size = mp.Vector3(Sx, Sy, Sz)\n", + "\n", + "substrate_thickness = Sz/2\n", + "substrate_size = mp.Vector3(design_region_width, design_region_length, substrate_thickness)\n", + "substrate_center = mp.Vector3(\n", + " 0,\n", + " 0,\n", + " - 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 = [\n", + " substrate_block,\n", + " design_region_block,\n", + "]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "607d05b7", + "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)]\n", + "#src = mp.ContinuousSource(frequency=fcen, fwidth=fwidth, is_integrated=True)\n", + "#source = [mp.Source(src, component=mp.Ex, size=source_size, center=source_center)]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cdbe5c87", + "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", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5f775457", + "metadata": {}, + "outputs": [], + "source": [ + "######MONITOR######\n", + "far_field_obs_point_z_coord = 15\n", + "far_z = [mp.Vector3(0,0,far_field_obs_point_z_coord)] #this is the far-field-focal-point\n", + "NearRegions = [\n", + " 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,Ey,Ez,Hx,Hy,Hz\n", + "FarFields = mpa.Near2FarFields(sim,NearRegions,far_z)\n", + "ob_list = [FarFields]\n", + "\n", + "#this the objective to function\n", + "#the objective is to maximize the obj function,\n", + "def J1(FF):\n", + " #in this case we want to maximize field intensity\n", + " # at the position of the focal spot \n", + " obj_to_return = npa.mean(npa.abs(FF[0,:,1]) ** 2)\n", + " return obj_to_return\n", + "#returns the first point of all frequencies of the Ey component" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "27fb5c57", + "metadata": {}, + "outputs": [], + "source": [ + "######CHECK PLOT######\n", + "output_plane = mp.Volume(\n", + " center=mp.Vector3(0,0,0), \n", + " size=mp.Vector3(Sx-2*pml_size,0,Sz-2*pml_size))\n", + "\n", + "output_plane_xy = mp.Volume(\n", + " center=mp.Vector3(0,0,design_region_height / 2),\n", + " size=mp.Vector3(Sx-2*pml_size,Sy-2*pml_size,0)\n", + ")\n", + "\n", + "fig = sim.plot2D(show_sources=True,show_monitors=True,show_epsilon=True, output_plane=output_plane)\n", + "#sim.run(until=10)\n", + "#sim.plot2D(fields=mp.Ex,output_plane=output_plane)\n", + "#plt.savefig(\"fields.pdf\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9c576fdb", + "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=500, #maximum optimizer run time\n", + ")\n", + "\n", + "######CHECK PLOT######\n", + "figure = opt.plot2D(True, output_plane=output_plane)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b315f21d", + "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) = plt.subplots(1, 2, figsize=(12, 5))\n", + " \n", + " opt.plot2D(\n", + " False, \n", + " ax=ax1,\n", + " output_plane=output_plane,\n", + " show_sources = True,\n", + " show_monitors = True,\n", + " show_boundary_layers = True,\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 = True,\n", + " show_monitors = True,\n", + " show_boundary_layers = True,\n", + " )\n", + " circ2 = Circle((2,2), minimum_length/2)\n", + " ax2.add_patch(circ2)\n", + " \n", + " if not os.path.exists(\"results\"):\n", + " os.makedirs(\"results\")\n", + "\n", + " plt.savefig(\"results/iter_{}.pdf\".format(cur_iter[0]+1))\n", + "\n", + " cur_iter[0] += 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_max_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": "164e57d7", + "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 = True,\n", + " show_monitors = True,\n", + " show_boundary_layers = True,\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 = True,\n", + " show_monitors = True,\n", + " show_boundary_layers = True,\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.pdf\")\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 +} diff --git a/mode_converter.py b/mode_converter.py new file mode 100644 index 0000000..69232bf --- /dev/null +++ b/mode_converter.py @@ -0,0 +1,355 @@ +# This is a first attempt at creating an inverse designed fiber mode converter in meep. At first, we try to implement and combine tutorials from both tidy3d and meep. + + +import meep as mp +import meep.adjoint as mpa +import numpy as np +from autograd import numpy as npa +from autograd import tensor_jacobian_product, grad +import nlopt +from matplotlib import pyplot as plt +from matplotlib.patches import Circle +import scipy as sp +import os + + +######PHYSICS SETUP###### + +mp.verbosity(1) + +subs_index = 3.4 +subs_perm = subs_index**2 + +Si = mp.Medium(index = subs_index) +Air = mp.Medium(index = 1.0) + +wavelength = 1.55 +nf = 3 #number of frequencies to test +frequencies = np.array([1/1.5, 1 / 1.55, 1 / 1.6]) + +#distance between PML and source/monitor +buffer = 2.5 + +######PML###### +pml_size = 1.0 +pml_layers = [mp.PML(pml_size)] + + +######FILTER###### +minimum_length = 0.09 # minimum length scale (microns) +eta_i = ( + 0.5 # blueprint (or intermediate) design field thresholding point (between 0 and 1) +) +eta_e = 0.55 # 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) + + +######DESIGN REGION###### +design_region_width = 4.0 +design_region_length = 4.0 +design_region_height = 1.0 + +resolution = 15 # pixels/um +design_region_resolution = int(resolution) + +#number of pixels in the design region +Nx = int(design_region_width * design_region_resolution) + 1 +Ny = int(design_region_length * design_region_resolution) + 1 + +design_variables = mp.MaterialGrid(mp.Vector3(Nx, Ny), Air, Si, grid_type="U_MEAN") + +design_region = mpa.DesignRegion( + design_variables, + volume=mp.Volume( + center = mp.Vector3(0,0,design_region_height / 2), + size = mp.Vector3(design_region_width, + design_region_length, + design_region_height), + ), +) + +######CELL GEOMETRY###### +Sx = 2 * pml_size + design_region_width +Sy = 2 * pml_size + design_region_length +Sz = 2 * pml_size + design_region_height + 2*buffer +print("Cell size: ", Sx, Sy, Sz) +cell_size = mp.Vector3(Sx, Sy, Sz) + +substrate_thickness = Sz/2 +substrate_size = mp.Vector3(design_region_width, design_region_length, substrate_thickness) +substrate_center = mp.Vector3( + 0, + 0, + - substrate_thickness / 2) + +substrate_block = mp.Block( + center = substrate_center, + size = substrate_size, + material = Si, +) + +design_region_block = mp.Block( + center = design_region.center, + size = design_region.size, + material = design_variables, +) + +geometry = [ + substrate_block, + design_region_block, +] + + +######SOURCE###### +fcen = 1 / wavelength +width = 0.2 +fwidth = width * fcen +source_center = [0, 0, -3] +source_size = mp.Vector3(design_region_width, design_region_length, 0) +src = mp.GaussianSource(frequency=fcen, fwidth=fwidth) +source = [mp.Source(src, component=mp.Ex, size=source_size, center=source_center)] +#src = mp.ContinuousSource(frequency=fcen, fwidth=fwidth, is_integrated=True) +#source = [mp.Source(src, component=mp.Ex, size=source_size, center=source_center)] + + +######SIMULATION###### +kpoint = mp.Vector3() +sim = mp.Simulation( + cell_size = cell_size, + boundary_layers = pml_layers, + geometry = geometry, + sources = source, + default_material = Air, + resolution = resolution, +) + + +######MONITOR###### +far_field_obs_point_z_coord = 15 +far_z = [mp.Vector3(0,0,far_field_obs_point_z_coord)] #this is the far-field-focal-point +NearRegions = [ + mp.Near2FarRegion( + center=mp.Vector3(0,0, design_region_height + wavelength), + size=mp.Vector3(design_region_width,design_region_length,0), + weight=+1, + ) +] +#mpa.Near2FarFields return a numpy array with shape (num_of_points,nfreq,6) where the +#third axis are the field components Ex,Ey,Ez,Hx,Hy,Hz +FarFields = mpa.Near2FarFields(sim,NearRegions,far_z) +ob_list = [FarFields] + +#this the objective to function +#the objective is to maximize the obj function, +def J1(FF): + #in this case we want to maximize field intensity + # at the position of the focal spot + obj_to_return = npa.mean(npa.abs(FF[0,:,1]) ** 2) + return obj_to_return +#returns the first point of all frequencies of the Ey component + + +######CHECK PLOT###### +output_plane = mp.Volume( + center=mp.Vector3(0,0,0), + size=mp.Vector3(Sx-2*pml_size,0,Sz-2*pml_size)) + +output_plane_xy = mp.Volume( + center=mp.Vector3(0,0,design_region_height / 2), + size=mp.Vector3(Sx-2*pml_size,Sy-2*pml_size,0) +) + +fig = sim.plot2D(show_sources=True,show_monitors=True,show_epsilon=True, output_plane=output_plane) +#sim.run(until=10) +#sim.plot2D(fields=mp.Ex,output_plane=output_plane) +#plt.savefig("fields.pdf") + + +######OPTIMIZATION PROBLEM###### +''' +An OptimizationProblem class object knows how to do one basic thing: +Given an input vector of design variables, compute the objective function value +(forward calculation) and optionally its gradient (adjoint calculation) +The gradient of the objective function is computed against its arguments, which are +the FarFields. +''' +opt = mpa.OptimizationProblem( + simulation=sim, #the optimizer uses this simulation object + objective_functions=[J1], #the optimizer call this objective function + objective_arguments=ob_list, #the optimizer passes this list of arguments to the objective function + design_regions=[design_region], #the optimizer updates the design variables in this design region + frequencies=frequencies, #the optimizer runs the simulation at these frequencies + maximum_run_time=500, #maximum optimizer run time +) + +######CHECK PLOT###### +figure = opt.plot2D(True, output_plane=output_plane) +plt.savefig("results/simulation_setup.pdf") + + + +def mapping(x,eta,beta): + ''' + Arguments: + x: design variables (grey scale) + eta: thresholding point for projection (between 0 and 1) + beta: binarization strength for projection (positive) + Returns: + projected_field: the filtered and projected design variables + ''' + #mapping is necessary to reach a physical design. + #we start from the design variables (grey scale), which are values in the range (0,1). + #typically, we set an initial value of 0.5. + #filtering makes a local point behave more like its neighbourhood, avoiding checker patterns + filtered_field = mpa.conic_filter( + x, + filter_radius, + design_region_width, + design_region_length, + [design_region_resolution], + ) + + #projection + #this is an hyperbolic tangent projection + #after filtering, the design is still blurry, so we project it + #using an tanh projection method to push the values towards the upper/lower bound + #If a filtered value is >0.5, then we push it to 1, otherwise, if <0.5, we push it to 0. + #parameter beta represents the degree of binarization, + #parameter eta is threshold point + projected_field = mpa.tanh_projection(filtered_field,beta,eta) + + #the following lines introduce symmetry. Need to understand if necessary + #projected_field = ( + # npa.flipud(projected_field) + projected_field + #)/2 + + #interpolate to actual materials + return projected_field.flatten() + +evaluation_history = [] +cur_iter = [0] + +def f(v, gradient, cur_beta): + print("current iteration: {}".format(cur_iter[0]+1)) + + #value and grad + f0, dJ_du = opt([mapping(v,eta_i,cur_beta)]) + if gradient.size > 0: + # NOTE: gradient is assigned in-place! No need to return. + gradient[...] = tensor_jacobian_product(mapping, 0)( + v, eta_i, cur_beta, np.sum(dJ_du, axis=1) + ) + + evaluation_history.append(np.real(f0)) + + fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5)) + + opt.plot2D( + False, + ax=ax1, + output_plane=output_plane, + show_sources = True, + show_monitors = True, + show_boundary_layers = True, + ) + circ1 = Circle((2,2), minimum_length/2) + ax1.add_patch(circ1) + + opt.plot2D( + False, + ax=ax2, + output_plane=output_plane_xy, + show_sources = True, + show_monitors = True, + show_boundary_layers = True, + ) + circ2 = Circle((2,2), minimum_length/2) + ax2.add_patch(circ2) + + if not os.path.exists("results"): + os.makedirs("results") + + plt.savefig("results/iter_{}.pdf".format(cur_iter[0]+1)) + + cur_iter[0] += cur_iter[0] + 1 + + return np.real(f0) + + +algorithm = nlopt.LD_MMA +n = Nx * Ny #number of parameters to optimize + +x = np.ones((n,)) * 0.5 #initial guess +#x = np.insert(x, 0, -1) +#lower and upper bounds for design variables +lb = np.zeros((Nx * Ny,)) +ub = np.ones((Nx * Ny,)) + +cur_beta = 4 #starting beta +beta_scale = 2 #beta scaling (update) factor +num_betas = 6 #number of betas to test +update_factor = 12 #number of iterations per beta +#total number of simulations = num_betas * update_factor * 2 (forward + adjoint) +ftol = 1e-5 + +optimization_start = True +if optimization_start: + for iters in range(num_betas): + solver = nlopt.opt(algorithm, n) + solver.set_lower_bounds(lb) + solver.set_upper_bounds(ub) + solver.set_max_objective(lambda v, g: f(v, g, cur_beta)) + solver.set_maxeval(update_factor) #number of iterations for each beta + solver.set_ftol_rel(ftol) + #observe that x is updated in place, so at + # the end it will contain the optimized parameters + x[...] = solver.optimize(x) + cur_beta = cur_beta * beta_scale #update beta + print("current beta: ", cur_beta) + + +######RESULTS###### +#plot the figure of merit over iterations +plt.figure() +plt.plot(evaluation_history) +plt.xlabel("Iteration") +plt.ylabel("Figure of Merit") +plt.grid() +plt.savefig("results/figure_of_merit.pdf") + +#plot and save the final design +np.savetxt("results/final_design.txt", mapping(x, eta_i, cur_beta)) +opt.update_design([mapping(x, eta_i, cur_beta)]) + +fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5)) + +opt.plot2D( + False, + ax=ax1, + output_plane=output_plane, + show_sources = True, + show_monitors = True, + show_boundary_layers = True, +) + +circ1 = Circle((2,2), minimum_length/2) +ax1.add_patch(circ1) +ax1.axis("off") + +opt.plot2D( + False, + ax=ax2, + output_plane=output_plane_xy, + show_sources = True, + show_monitors = True, + show_boundary_layers = True, +) +circ2 = Circle((2,2), minimum_length/2) +ax2.add_patch(circ2) +ax2.axis("off") + +plt.savefig("results/final_design.pdf") + + +