-
Notifications
You must be signed in to change notification settings - Fork 76
/
Copy pathoptimize_1_3.py
185 lines (156 loc) · 7.13 KB
/
optimize_1_3.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
import argparse
import numpy as np
import autograd.numpy as npa
import matplotlib.pylab as plt
from autograd.scipy.signal import convolve as conv
from skimage.draw import circle
import ceviche
from ceviche import fdfd_ez, jacobian
from ceviche.optimizers import adam_optimize
from ceviche.utils import imarr, get_value
from ceviche.modes import insert_mode
import collections
# Create a container for our slice coords to be used for sources and probes
Slice = collections.namedtuple('Slice', 'x y')
def init_domain(Nx, Ny, Npml, space=10, wg_width=10, space_slice=5):
"""Initializes the domain and design region
space : The space between the PML and the structure
wg_width : The feed and probe waveguide width
space_slice : The added space for the probe and source slices
"""
rho = np.zeros((Nx, Ny))
design_region = np.zeros((Nx, Ny))
# Input probe slice
input_slice = Slice(x=np.array(Npml+1), y=np.arange(int(Ny/2-wg_width/2-space_slice), int(Ny/2+wg_width/2+space_slice)))
output_slices = []
# Input waveguide
rho[0:int(Npml+space),int(Ny/2-wg_width/2):int(Ny/2+wg_width/2)] = 1
# Output waveguide 1
rho[int(Nx-Npml-space)::,int(Npml+space):int(Npml+space+wg_width)] = 1
output_slices.append(Slice(x=np.array(Nx-Npml-1), y=np.arange(int(Npml+space-space_slice), int(Npml+space+wg_width+space_slice))))
# Output waveguide 2
rho[int(Nx-Npml-space)::,int(Ny/2-wg_width/2):int(Ny/2+wg_width/2)] = 1
output_slices.append(Slice(x=np.array(Nx-Npml-1), y=np.arange(int(Ny/2-wg_width/2-space_slice), int(Ny/2+wg_width/2+space_slice))))
# Output waveguide 3
rho[int(Nx-Npml-space)::,int(Ny-Npml-space-wg_width):int(Ny-Npml-space)] = 1
output_slices.append(Slice(x=np.array(Nx-Npml-1), y=np.arange(int(Ny-Npml-space-wg_width-space_slice), int(Ny-Npml-space+space_slice))))
design_region[Npml+space:Nx-Npml-space, Npml+space:Ny-Npml-space] = 1
rho[Npml+space:Nx-Npml-space, Npml+space:Ny-Npml-space] = 0.5
return rho, design_region, input_slice, output_slices
def operator_proj(rho, eta=0.5, beta=100):
"""Density projection
"""
return npa.divide(npa.tanh(beta * eta) + npa.tanh(beta * (rho - eta)), npa.tanh(beta * eta) + npa.tanh(beta * (1 - eta)))
def operator_blur(rho, radius=2):
"""Blur operator implemented via two-dimensional convolution
"""
rr, cc = circle(radius, radius, radius+1)
kernel = np.zeros((2*radius+1, 2*radius+1), dtype=np.float)
kernel[rr, cc] = 1
kernel=kernel/kernel.sum()
# For whatever reason HIPS autograd doesn't support 'same' mode, so we need to manually crop the output
return conv(rho, kernel, mode='full')[radius:-radius,radius:-radius]
def make_rho(rho, design_region, radius=2):
"""Helper function for applying the blur to only the design region
"""
lpf_rho = operator_blur(rho, radius=radius) * design_region
bg_rho = rho * (design_region==0).astype(np.float)
return bg_rho + lpf_rho
def viz_sim(epsr):
"""Solve and visualize a simulation with permittivity 'epsr'
"""
simulation = fdfd_ez(omega, dl, epsr, [Npml, Npml])
Hx, Hy, Ez = simulation.solve(source)
fig, ax = plt.subplots(1, 2, constrained_layout=True, figsize=(6,3))
ceviche.viz.real(Ez, outline=epsr, ax=ax[0], cbar=False)
ax[0].plot(input_slice.x*np.ones(len(input_slice.y)), input_slice.y, 'g-')
for output_slice in output_slices:
ax[0].plot(output_slice.x*np.ones(len(output_slice.y)), output_slice.y, 'r-')
ceviche.viz.abs(epsr, ax=ax[1], cmap='Greys');
plt.show()
return (simulation, ax)
""" DEFINE ALL THE PARAMETERS """
parser = argparse.ArgumentParser()
# Number of epochs in the optimization
parser.add_argument('-Nsteps', default=50, type=int)
# Step size for the Adam optimizer
parser.add_argument('-step_size', default=1e0, type=float)
# Angular frequency of the source in 1/s
parser.add_argument('-omega', default=2*np.pi*200e12, type=float)
# Spatial resolution in meters
parser.add_argument('-dl', default=50e-9, type=float)
# Number of pixels in x-direction
parser.add_argument('-Nx', default=110, type=int)
# Number of pixels in y-direction
parser.add_argument('-Ny', default=130, type=int)
# Number of pixels in the PMLs in each direction
parser.add_argument('-Npml', default=20, type=int)
# Minimum value of the relative permittivity
parser.add_argument('-epsr_min', default=1.0, type=float)
# Maximum value of the relative permittivity
parser.add_argument('-epsr_max', default=12.0, type=float)
# Radius of the smoothening features
parser.add_argument('-blur_radius', default=3, type=int)
# Strength of the binarizing projection
parser.add_argument('-beta', default=100.0, type=float)
# Middle point of the binarizing projection
parser.add_argument('-eta', default=0.5, type=float)
# Space between the PMLs and the design region (in pixels)
parser.add_argument('-space', default=10, type=int)
# Width of the waveguide (in pixels)
parser.add_argument('-wg_width', default=10, type=int)
# Length in pixels of the source/probe slices on each side of the center point
parser.add_argument('-space_slice', default=7, type=int)
# Port to which transmission will be maximized
parser.add_argument('-objective_port_ind', default=0, type=int)
args = parser.parse_args()
Nsteps = args.Nsteps
step_size = args.step_size
omega = args.omega
dl = args.dl
Nx = args.Nx
Ny = args.Ny
Npml = args.Npml
epsr_min = args.epsr_min
epsr_max = args.epsr_max
blur_radius = args.blur_radius
beta = args.beta
eta = args.eta
space = args.space
wg_width = args.wg_width
space_slice = args.space_slice
objective_port_ind = args.objective_port_ind
""" END PARAMETER DEFINITION """
# Setup initial structure
rho_init, design_region, input_slice, output_slices = \
init_domain(Nx, Ny, Npml, space=space, wg_width=wg_width, space_slice=space_slice)
epsr = epsr_min + (epsr_max-epsr_min) * make_rho(rho_init, design_region, radius=blur_radius)
# Setup source
source = insert_mode(omega, dl, input_slice.x, input_slice.y, epsr)
# Setup probes
probes = []
for output_slice in output_slices:
probe = insert_mode(omega, dl, output_slice.x, output_slice.y, epsr)
probes.append(probe)
# Simulate initial device
(simulation, ax) = viz_sim(epsr)
# Define optimization objective
def measure_modes(Ez, probe_ind=0):
return npa.abs(npa.sum(npa.conj(Ez)*probes[probe_ind]))
def objective(rho, probe_ind=objective_port_ind):
rho = rho.reshape((Nx, Ny))
_rho = make_rho(rho, design_region, radius=blur_radius)
epsr = epsr_min + (epsr_max-epsr_min)*operator_proj(_rho, beta=beta, eta=eta) * design_region
simulation.eps_r = epsr
_, _, Ez = simulation.solve(source)
return measure_modes(Ez, probe_ind=probe_ind)
# Run optimization
objective_jac = jacobian(objective, mode='reverse')
(rho_optimum, loss) = adam_optimize(objective, rho_init.flatten(), objective_jac,
Nsteps=Nsteps, direction='max',
step_size=step_size)
rho_optimum = rho_optimum.reshape((Nx, Ny))
# Simulate optimal device
epsr = epsr_min + (epsr_max-epsr_min)*operator_proj(make_rho(rho_optimum,
design_region, radius=blur_radius), beta=beta, eta=eta)
(simulation, ax) = viz_sim(epsr)