-
Notifications
You must be signed in to change notification settings - Fork 3
/
1d-theano.py
executable file
·169 lines (136 loc) · 5.04 KB
/
1d-theano.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
#!/usr/bin/env python3
# 1D FDTD, Ey/Hx mode
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
from matplotlib import animation
from scipy import signal
import time
import theano
import theano.sparse
import theano.tensor as T
# TODO:
# * Receive hints and simulation parameters as command-line parameters
# Basic constants
m0 = 4*np.pi*10**-7 # N/A**2
e0 = 8.854187817*10**-12 # F/m
c0 = 1.0/np.sqrt(m0*e0) # The speed of light
# Simulation parameters
KHz = 1000.0
MHz = 1000000.0
GHz = 1000000000.0
ps = 10**-12
ns = 10**-9
us = 10**-6
# Material layers
# Layers are overwriting each other, last write wins
# First layer should be free space.
# (mr, er, start, width)
layers = [
(1.0, 1.0, 0.0, 1.0), # Free space
(1.0, 2.5, 0.6, 0.2) # A 20cm slab of plastic
]
# Calculate maximal refractive index
n_max = 1.0
n_min = 1000000.0 # TODO: init to first layer
for layer in layers:
n = layer[0] / layer[1]
if n > n_max:
n_max = n
if n < n_min:
n_min = n
space_size = 1.0 # meters
freq_max = 100*GHz # maximal resolvable frequency
lamb_min = c0 / (freq_max * n_max) # minimal wavelength
dzpmwl = 10 # delta-z per minimal wavelength, a rule-of-thumb constant
dz = lamb_min / dzpmwl # Spatial step size, meters
gridsize = int(space_size / dz) # Size of the grid in cells
simlen = 5 * space_size / c0 # Simulation length, seconds (5 travels back & forth)
dt = n_min * dz / (2*c0) # From the Courant-Friedrichs-Lewy condition. This is a rule of thumb
steps = int(simlen / dt) # Number of simulation steps
print("simulation length:", simlen)
print("grid size:", gridsize)
print("steps:", steps)
print("dt:", dt)
print("dz:", dz)
mr = np.ones(gridsize, dtype=np.float32) # permeability, can be diagonally anisotropic
er = np.ones(gridsize, dtype=np.float32) # permittivity, can be diagonally anisotropic
for layer in layers:
# TODO: snap layers to grid / snap grid to layers?
for i in range(max(0, int(layer[2]/dz)), min(gridsize, int((layer[2]+layer[3])/dz))):
er[i] = layer[0]
mr[i] = layer[1]
# Update coeffitients, using normalized magnetic field
mkhx = c0*dt/mr
mkey = c0*dt/er
# Yee grid scheme
# dx, dy, dz, must be as square as possible, but can be different
# Function values are assigned at the middle of the square
#
# Field components are staggered at different places around
# the grid unit cube / square
#
# * This helps naturally satisfy the divergence equations
# * Material boundaries are naturally handled
# * Easier to calculate discrete curls
# * WARNING: field components can be in different materials!
E = theano.shared(np.zeros(gridsize, dtype=np.float32)) # Electric field
H = theano.shared(np.zeros(gridsize, dtype=np.float32)) # Normalized magnetic field
# Display
fig = plt.figure()
ax = plt.axes(ylim=(-5, 5))
line1, = ax.plot(np.linspace(0.0, space_size, gridsize), np.zeros(gridsize, dtype=np.float32), 'r-')
line2, = ax.plot(np.linspace(0.0, space_size, gridsize), np.zeros(gridsize, dtype=np.float32), 'b-')
for layer in layers:
n = layer[0]/layer[1]
ax.add_patch(patches.Rectangle((layer[2], -20), layer[3], 40, color=(n, n, n)))
# Sinc function source
def sinc_source(er, ur, period, t0, t):
a_corr = -np.sqrt(er/ur) # amplitude correction term
t_corr = np.sqrt(er*ur)*dz/(2*c0) + dt/2 # Time correction term
return (
# H field
a_corr * np.sinc((t-t0)*2/period + t_corr),
# E field
np.sinc((t-t0)*2/period) # E
)
# Gaussian pulse source
def gausspulse_source(er, ur, t0, tau, t):
a_corr = -np.sqrt(er/ur) # amplitude correction term
t_corr = np.sqrt(er*ur)*dz/(2*c0) + dt/2 # Time correction term
return (
a_corr * np.exp(-((t-t0)/tau)**2 + t_corr),
np.exp(-((t-t0)/tau)**2)
)
# Outputs 1.0 at time 0
def blip_source(t):
return (0.0, 1.0) if t == 0 else (0.0, 0.0)
def init_animation():
global line1, line2
return line1, line2
src = T.fscalar()
step1 = theano.function([src], None, updates=[
(H, T.inc_subtensor(T.inc_subtensor(H[1:-1], mkhx[1:-1] * (E[2:] - E[1:-1]) / dz)[500:501], [src]))
])
step2 = theano.function([src], None, updates=[
(E, T.inc_subtensor(T.inc_subtensor(E[1:-1], mkey[1:-1] * (H[1:-1] - H[:-2]) / dz)[500:501], [src]))
])
i = 0
def animate(_):
global i, ax, line1, line2, H, E, mkhx, mkey
print(i)
for i in range(i, i+100):
t = i*dt
src_ = gausspulse_source(1.0, 1.0, 200*ps, 50*ps, t)
step1(np.float32(src_[0]))
step2(np.float32(src_[1]))
# Simply dampen at the edges instead of using the messy perfect edge or PML method.
#E[:100] *= np.linspace(0.985,1.0,100)
#H[:100] *= np.linspace(0.985,1.0,100)
#E[-100:] *= np.linspace(1.0,0.985,100)
#H[-100:] *= np.linspace(1.0,0.985,100)
line1.set_ydata(E.get_value())
line2.set_ydata(H.get_value())
return line1, line2
anim = animation.FuncAnimation(fig, animate, init_func=init_animation, interval=0, blit=True)
plt.show()