developer.html
来自「麻省理工的计算光子晶体的程序」· HTML 代码 · 共 275 行
HTML
275 行
<!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 iscomputing, and then describe how the program and computation arebroken up into different portions of the code.<p>Forgive the primitive math typography below; this will be rectifiedwhen <a href="http://www.w3.org/Math/">MathML</a> is supported in a <ahref="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 ofphotonic band structure calculations and the algorithms that weemploy. 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>PhotonicCrystals: 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/abstract.cfm?URI=OPEX-8-3-173">Block-iterativefrequency-domain methods for Maxwell's equations in a planewavebasis</a>," <i>Optics Express</i> <b>8</b>, no. 3, 173-190 (2001).</ul><p>The MIT Photonic-Bands Package takes a periodic dielectricstructure and computes the <em>eigenmodes</em> of that structure,which are the electromagnetic waves that can propagate through thestructure with a definite frequency. This corresponds to solving aneigenvalue 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 theMaxwell operator <code>curl 1/epsilon curl</code>. We also have anadditional constraint, that <code>div h</code> be zero (the magneticfield must be "transverse").<p>Since the structure is periodic, we can also invoke Bloch's theoremto write the states in the form <code>exp(i k*x)</code> times aperiodic function, where <code>k</code> is the Bloch wavevector. So,at each k-point (Bloch wavevector), we need to solve for a discreteset of eigenstates, the photonic bands of the structure.<p>To solve for the eigenstates on a computer, we must expand themagnetic field in some basis, where we truncate the basis to somefinite number of points to discretize the problem. For example, wecould use a traditional finite-element basis in which the field istaken on a finite number of mesh points and linearly interpolated inbetween. However, it is expensive to enforce the transversalityconstraint in this basis. Instead, we use a Fourier (spectral) basis,expanding the periodic part of the field in terms of <code>exp(iG*x)</code> planewaves. In this basis, the transversality constraintis easy to maintain, as it merely implies that the planewaveamplitudes must be orthogonal to <code>k + G</code>.<p>In order to find the eigenfunctions, we could compute the elementsof <code>M</code> explicitly in our basis, and then call LAPACK orsome similar code to find the eigenvectors and eigenvalues. For athree-dimensional calculation, this could mean finding theeigenvectors of a matrix with hundreds of thousands of elements on aside--daunting merely to store, much less compute. Fortunately, weonly want to know a few eigenvectors, not hundreds of thousands, so wecan use much less expensive <em>iterative</em> methods that don'trequire us to store <code>M</code> explicitly.<p>Iterative eigensolvers require only that one supply a routine tooperate <code>M</code> on a vector (function). Starting with aninitial guess for the eigenvector, they then converge quickly to theactual eigenvector, stopping when the desired tolerance is achieved.There are many iterative eigensolver methods; we use a preconditionedblock minimization of the Rayleigh quotient which is further describedin the file <code>src/matrices/eigensolver.c</code>. In the Fourierbasis, applying <code>M</code> to a function is relatively easy: thecurls become cross products with <code>i (k + G)</code>; themultiplication by <code>1/epsilon</code> is performed by using an FFTto transform to the spatial domain, multiplying, and then transformingback with an inverse FFT. For more information and references oniterative eigensolvers, see the paper cited above.<p>We also support a "targeted" eigensolver. A typical iterativeeigensolver finds the <code>p</code> lowest eigenvalues andeigenvectors. Instead, we can find the <code>p</code> eigenvaluesclosest to a given frequency <code>w0</code> by solving for theeigenvalues 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 thoseclosest to <code>w0</code> the smallest.<p>The eigensolver we use is preconditioned, which means that convergencecan be greatly improved by suppling a good preconditioner matrix.Finding a good preconditioner involves making an approximate inverseof <code>M</code>, and is something of a black art with lots of trialand error.<h2><a name="epsilon">Dielectric Function Computation</a></h2><p>The initialization of the dielectric function deserves someadditional discussion, both because it is crucial for goodconvergence, and because we use somewhat complicated algorithms forperformance reasons.<p>To ameliorate the convergence problems caused in a planewave basisby a discontinuous dielectric function, the dielectric function issmoothed (averaged) at the resolution of the grid. Another way ofthinking about it is that this brings the average dielectric constant(over the grid) closer to its true value. Since differentpolarizations of the field prefer different averaging methods, one hasto construct an effective dielectric tensor at the boundaries betweendielectrics, as described by the paper referenced above.<p>This averaging has two components. First, at each grid point thedielectric constant (epsilon) and its inverse are averaged over auniform mesh extending halfway to the neighboring grid points. (Themesh resolution is controlled by the <code>mesh-size</code> user inputvariable.) Second, for grid points on the boundary between twodielectrics, we compute the vector normal to the dielectric interface;this is done by averaging the "dipole moment" of the dielectricfunction over a spherically-symmetric distribution of points. Thenormal vector and the two averages of epsilon are then combined intoan 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 functionepsilon(<b>r</b>), which returns the dielectric constant for a givenposition <b>r</b>. This epsilon function must be as efficient aspossible, because it is evaluated a large number of times: the size ofthe grid multiplied by <code>mesh-size</code><sup>3</sup> (in threedimensions).<p>To specify the geometry, the user provides a list of geometricobjects (blocks, spheres, cylinders and so on). These are parsed intoan efficient data structure and are used to to provide the epsilonfunction described above. (All of this is handled by the libctlgeomcomponent of libctl, described below.) At the heart of the epsilonfunction is a routine to return the geometric object enclosing a givenpoint, taking into account the fact that the objects are periodic inthe lattice vectors. Our first algorithm for doing this was a simplelinear search through the list of objects and their translations bythe lattice vectors, but this proved to be too slow, especially insupercell calculations where there are many objects. We addressed theperformance problem in two ways. First, for each object we constructa bounding box, with which point inclusion can be tested rapidly.Second, we build a hierarchical tree of bounding boxes, recursivelypartitioning the set of objects in the cell. This allows us to searchfor the object containing a point in a time logarithmic in the numberof 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 ofthe user interface, and to keep the eigensolver routines independentof the operator they are computing the eigenvector of. Thecomputational code is located in the <code>src/</code> directory, witha few major subdirectories, described below. The Guile-based userinterface 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 itreturns 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 moreinformation about the data structures. In particular, you should usethe supplied functions (<code>create_evectmatrix</code>, etcetera) tocreate 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 andsingle 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 supportboth real and complex numbers elsewhere in the code.<p>Also in this directory is <code>blasglue.c</code>, a set of wrapperroutines to make it convienient to call BLAS and LAPACK routines fromC instead of Fortran.<h3>src/util/</h3><p>As its name implies, this is simply a number of utility routinesfor 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 thecode to improve robustness. There are also debugging versions ofmalloc/free (which perform lots of paranoia tests, enabled by<code>--enable-debug-malloc</code> in <code>configure</code>), and MPIglue routines that allow the program to operate without the MPIlibraries.<h3>src/matrixio</h3><p>This section contains code to abstract I/O for eigenvectors andsimilar matrices, providing a simpler layer on top of the HDF5interface. This could be modified to support HDF4 or other I/Oformats.<h3>src/maxwell/</h3><p>The <code>maxwell/</code> directory contains all knowledge ofMaxwell's equations used by the program. It implements functions toapply 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 atransverse 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>. Onemust also initialize the dielectric function using<code>set_maxwell_dielectric</code> by supplying a function thatreturns the dielectric constant for any given coordinate. You canalso restrict yourself to TE or TM polarizations in two dimensions bycalling <code>set_maxwell_data_polarization</code>.<p>This directory also contains functions<code>maxwell_compute_dfield</code>, etcetera, to compute theposition-space fields from the Fourier-transform representationreturned 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. Thismeans that the user-interface code (in <code>mpb.c</code>) is fairlyshort, consisting of a number of small functions that are callable bythe 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 <ahref="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 parameterslike 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 subroutinesthat are visible to the Guile user. It also defines a number ofScheme subroutines for the user to call directly, like<code>(run)</code>. It is often simpler and more flexible to definefunctions like this in Scheme rather than in C.<p>All of the code to handle the geometric objects resides inlibctlgeom, a set of Scheme and C utility functions included withlibctl (see the file <code>utils/README</code> in the libctl package).(These functions could also be useful in other programs, such as atime-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>⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?