mirror of
https://github.com/edoverg/modeconverter.git
synced 2026-10-14 23:06:39 +00:00
first commit
This commit is contained in:
edoardo
@@ -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
|
||||
}
|
||||
@@ -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")
|
||||
|
||||
|
||||
|
||||
Reference in New Issue
Block a user