-
Notifications
You must be signed in to change notification settings - Fork 83
/
developer.html
274 lines (230 loc) · 13.4 KB
/
developer.html
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
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 3.2//EN">
<HTML><HEAD>
<TITLE>Developer Information</TITLE>
<LINK rel="Contents" href="index.html">
<LINK rel="Copyright" href="license.html">
<LINK rel="Start" href="index.html">
</HEAD>
<BODY TEXT="#000000" BGCOLOR="#FFFFFF">
Go to the <a href="acknowledgments.html">next</a>, <a href="user-ref.html">previous</a>, or <a href="index.html">main</a> section.
<hr>
<h1>Developer Information</h1>
<p>Here, we begin with a brief overview of what the program is
computing, and then describe how the program and computation are
broken up into different portions of the code.
<p>Forgive the primitive math typography below; this will be rectified
when <a href="http://www.w3.org/Math/">MathML</a> is supported in a <a
href="http://www.mozilla.org/projects/mathml/">decent browser</a>.
<h2><a name="math">The Mathematics of MPB</a></h2>
<p>This section provides a whirlwind tour of the mathematics of
photonic band structure calculations and the algorithms that we
employ. For more detailed information, see:
<ul>
<li><a title="the book at BarnesandNoble.com"
href="http://shop.barnesandnoble.com/booksearch/isbnInquiry.asp?isbn=0691037442"><i>Photonic
Crystals: Molding the Flow of Light</i></a>, by J. D. Joannopoulos,
R. D. Meade, and J. N. Winn (Princeton, 1995).
<li>Steven G. Johnson and J. D. Joannopoulos, "<a href="http://www.opticsexpress.org/oearchive/source/27937.htm">Block-iterative
frequency-domain methods for Maxwell's equations in a planewave
basis</a>," <i>Optics Express</i> <b>8</b>, no. 3, 173-190 (2001).
</ul>
<p>The MIT Photonic-Bands Package takes a periodic dielectric
structure and computes the <em>eigenmodes</em> of that structure,
which are the electromagnetic waves that can propagate through the
structure with a definite frequency. This corresponds to solving an
eigenvalue problem <code>M h = (w/c)^2 h</code>, where <code>h</code>
is the magnetic field, <code>w</code> is the frequency, and M is the
Maxwell operator <code>curl 1/epsilon curl</code>. We also have an
additional constraint, that <code>div h</code> be zero (the magnetic
field must be "transverse").
<p>Since the structure is periodic, we can also invoke Bloch's theorem
to write the states in the form <code>exp(i k*x)</code> times a
periodic function, where <code>k</code> is the Bloch wavevector. So,
at each k-point (Bloch wavevector), we need to solve for a discrete
set of eigenstates, the photonic bands of the structure.
<p>To solve for the eigenstates on a computer, we must expand the
magnetic field in some basis, where we truncate the basis to some
finite number of points to discretize the problem. For example, we
could use a traditional finite-element basis in which the field is
taken on a finite number of mesh points and linearly interpolated in
between. However, it is expensive to enforce the transversality
constraint in this basis. Instead, we use a Fourier (spectral) basis,
expanding the periodic part of the field in terms of <code>exp(i
G*x)</code> planewaves. In this basis, the transversality constraint
is easy to maintain, as it merely implies that the planewave
amplitudes must be orthogonal to <code>k + G</code>.
<p>In order to find the eigenfunctions, we could compute the elements
of <code>M</code> explicitly in our basis, and then call LAPACK or
some similar code to find the eigenvectors and eigenvalues. For a
three-dimensional calculation, this could mean finding the
eigenvectors of a matrix with hundreds of thousands of elements on a
side--daunting merely to store, much less compute. Fortunately, we
only want to know a few eigenvectors, not hundreds of thousands, so we
can use much less expensive <em>iterative</em> methods that don't
require us to store <code>M</code> explicitly.
<p>Iterative eigensolvers require only that one supply a routine to
operate <code>M</code> on a vector (function). Starting with an
initial guess for the eigenvector, they then converge quickly to the
actual eigenvector, stopping when the desired tolerance is achieved.
There are many iterative eigensolver methods; we use a preconditioned
block minimization of the Rayleigh quotient which is further described
in the file <code>src/matrices/eigensolver.c</code>. In the Fourier
basis, applying <code>M</code> to a function is relatively easy: the
curls become cross products with <code>i (k + G)</code>; the
multiplication by <code>1/epsilon</code> is performed by using an FFT
to transform to the spatial domain, multiplying, and then transforming
back with an inverse FFT. For more information and references on
iterative eigensolvers, see the paper cited above.
<p>We also support a "targeted" eigensolver. A typical iterative
eigensolver finds the <code>p</code> lowest eigenvalues and
eigenvectors. Instead, we can find the <code>p</code> eigenvalues
closest to a given frequency <code>w0</code> by solving for the
eigenvalues of <code>(M - (w0/c)^2)^2</code> instead of
<code>M</code>. This new operator has the same eigenvectors as
<code>M</code>, but its eigenvalues have been shifted to make those
closest to <code>w0</code> the smallest.
<p>The eigensolver we use is preconditioned, which means that convergence
can be greatly improved by suppling a good preconditioner matrix.
Finding a good preconditioner involves making an approximate inverse
of <code>M</code>, and is something of a black art with lots of trial
and error.
<h2><a name="epsilon">Dielectric Function Computation</a></h2>
<p>The initialization of the dielectric function deserves some
additional discussion, both because it is crucial for good
convergence, and because we use somewhat complicated algorithms for
performance reasons.
<p>To ameliorate the convergence problems caused in a planewave basis
by a discontinuous dielectric function, the dielectric function is
smoothed (averaged) at the resolution of the grid. Another way of
thinking about it is that this brings the average dielectric constant
(over the grid) closer to its true value. Since different
polarizations of the field prefer different averaging methods, one has
to construct an effective dielectric tensor at the boundaries between
dielectrics, as described by the paper referenced above.
<p>This averaging has two components. First, at each grid point the
dielectric constant (epsilon) and its inverse are averaged over a
uniform mesh extending halfway to the neighboring grid points. (The
mesh resolution is controlled by the <code>mesh-size</code> user input
variable.) Second, for grid points on the boundary between two
dielectrics, we compute the vector normal to the dielectric interface;
this is done by averaging the "dipole moment" of the dielectric
function over a spherically-symmetric distribution of points. The
normal vector and the two averages of epsilon are then combined into
an effective dielectric tensor for the grid point.
<p>All of this averaging is handled by a subroutine in
<code>src/maxwell/</code> (see below) that takes as input a function
epsilon(<b>r</b>), which returns the dielectric constant for a given
position <b>r</b>. This epsilon function must be as efficient as
possible, because it is evaluated a large number of times: the size of
the grid multiplied by <code>mesh-size</code><sup>3</sup> (in three
dimensions).
<p>To specify the geometry, the user provides a list of geometric
objects (blocks, spheres, cylinders and so on). These are parsed into
an efficient data structure and are used to to provide the epsilon
function described above. (All of this is handled by the libctlgeom
component of libctl, described below.) At the heart of the epsilon
function is a routine to return the geometric object enclosing a given
point, taking into account the fact that the objects are periodic in
the lattice vectors. Our first algorithm for doing this was a simple
linear search through the list of objects and their translations by
the lattice vectors, but this proved to be too slow, especially in
supercell calculations where there are many objects. We addressed the
performance problem in two ways. First, for each object we construct
a bounding box, with which point inclusion can be tested rapidly.
Second, we build a hierarchical tree of bounding boxes, recursively
partitioning the set of objects in the cell. This allows us to search
for the object containing a point in a time logarithmic in the number
of objects (instead of linear as before).
<h2><a name="code">Code Organization</a></h2>
<p>The code is organized to keep the core computation independent of
the user interface, and to keep the eigensolver routines independent
of the operator they are computing the eigenvector of. The
computational code is located in the <code>src/</code> directory, with
a few major subdirectories, described below. The Guile-based user
interface is completely contained within the <code>mpb-ctl/</code>
directory.
<h3>src/matrices/</h3>
<p>This directory contains the eigensolver, in
<code>eigensolver.c</code>, to which you pass an operator and it
returns the eigenvectors. Eigenvectors are stored using the
<code>evectmatrix</code> data structure, which holds <code>p</code>
eigenvectors of length <code>n</code>, potentially distributed over
<code>n</code> in MPI. See <code>src/matrices/README</code> for more
information about the data structures. In particular, you should use
the supplied functions (<code>create_evectmatrix</code>, etcetera) to
create and manipulate the data structures, where possible.
<p>The type of the eigenvector elements is determined by
<code>scalar.h</code>, which sets whether they are real or complex and
single or double precision. This is, in turn, controlled by the
<code>--disable-complex</code> and <code>--enable-single</code>
parameters to the <code>configure</code> script at install-time.
<code>scalar.h</code> contains macros to make it easier to support
both real and complex numbers elsewhere in the code.
<p>Also in this directory is <code>blasglue.c</code>, a set of wrapper
routines to make it convienient to call BLAS and LAPACK routines from
C instead of Fortran.
<h3>src/util/</h3>
<p>As its name implies, this is simply a number of utility routines
for use elsewhere in the code. Of particular note is
<code>check.h</code>, which defines a <code>CHECK(<i>condition</i>,
<i>error-message</i>)</code> macro that is used extensively in the
code to improve robustness. There are also debugging versions of
malloc/free (which perform lots of paranoia tests, enabled by
<code>--enable-debug-malloc</code> in <code>configure</code>), and MPI
glue routines that allow the program to operate without the MPI
libraries.
<h3>src/matrixio</h3>
<p>This section contains code to abstract I/O for eigenvectors and
similar matrices, providing a simpler layer on top of the HDF5
interface. This could be modified to support HDF4 or other I/O
formats.
<h3>src/maxwell/</h3>
<p>The <code>maxwell/</code> directory contains all knowledge of
Maxwell's equations used by the program. It implements functions to
apply the Maxwell operator to a vector (in <code>maxwell_op.c</code>)
and compute a good preconditioner (in <code>maxwell_pre.c</code>).
These functions operate upon a representation of the fields in a
transverse Fourier basis.
<p>In order to use these functions, one must first initialize a
<code>maxwell_data</code> structure with
<code>create_maxwell_data</code> (defined in <code>maxwell.c</code>)
and specify a k point with <code>update_maxwell_data_k</code>. One
must also initialize the dielectric function using
<code>set_maxwell_dielectric</code> by supplying a function that
returns the dielectric constant for any given coordinate. You can
also restrict yourself to TE or TM polarizations in two dimensions by
calling <code>set_maxwell_data_polarization</code>.
<p>This directory also contains functions
<code>maxwell_compute_dfield</code>, etcetera, to compute the
position-space fields from the Fourier-transform representation
returned by the eigensolver.
<h3>mpb-ctl/</h3>
<p>Here is the Guile-based user interface code for the eigensolver.
Instead of using Guile directly, this code is built on top of the
<code>libctl</code> library as described in previous sections. This
means that the user-interface code (in <code>mpb.c</code>) is fairly
short, consisting of a number of small functions that are callable by
the user from Guile.
<p>The core of the user interface is the file <code>mpb.scm</code>,
the <i>specifications file</i> for libctl as described in the <a
href="http://ab-initio.mit.edu/libctl/doc/">libctl manual</a>.
Actually, <code>mpb.scm</code> is generated by <code>configure</code>
from <code>mpb.scm.in</code> (in order to substitute in parameters
like the location of the libctl library); you should only edit
<code>mpb.scm.in</code> directly. (You can regenerate
<code>mpb.scm</code> simply by running <code>./config.status</code>
instead of re-running <code>configure</code>.)
<p>The specifications file defines the data structures and subroutines
that are visible to the Guile user. It also defines a number of
Scheme subroutines for the user to call directly, like
<code>(run)</code>. It is often simpler and more flexible to define
functions like this in Scheme rather than in C.
<p>All of the code to handle the geometric objects resides in
libctlgeom, a set of Scheme and C utility functions included with
libctl (see the file <code>utils/README</code> in the libctl package).
(These functions could also be useful in other programs, such as a
time-domain Maxwell's equation simulator.)
<hr>
Go to the <a href="acknowledgments.html">next</a>, <a href="user-ref.html">previous</a>, or <a href="index.html">main</a> section.
</BODY>
</HTML>