mirror of
https://github.com/edoverg/modeconverter.git
synced 2026-10-14 23:06:39 +00:00
464 lines
15 KiB
Plaintext
464 lines
15 KiB
Plaintext
{
|
|
"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",
|
|
"#setup a results folder\n",
|
|
"if not os.path.exists(\"results\"):\n",
|
|
" os.makedirs(\"results\")\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
|
|
}
|