CURLY3D: A finite element code for solving photonics problems in 3-D
Description
Curly3d solves the stationary Maxwell equations in the unit cell of a triply periodic, infinite domain containing dielectric and/or conducting material. The dielectric can have a non-positive definite permittivity tensor (epsilon) and/or a magnetic permeability tensor (mu). The finite element method is used to discretize the equation, yielding a generalized eigenvalue problem of the form A x = lambda B x. The basis functions are van Welij vector elements, with the unknown (scalar) amplitudes representing the field along the edges of straight hexahedral (brick) elements. To accelerate the numerical eigenvalue convergence, a mixed finite-element, finite-difference scheme is applied for the discretization of the identity operator.
Curly3d can solve either for the electric E or magnetic field H. The equation for H is similar to E except for the permittivity and permeability playing opposite roles. Note, however, that the E and H fields satisfy different boundary conditions at conductor interfaces.
A more extensive description of the numercal scheme used in Curly3d can be found in the proceedings of the ETOPIM 2002, Snowbird, Utah conference (July 15-19 2002).
Usage
Curly3d is a collection of Python modules (scripts). In order to use Curly3d, you'll need to write a driver, ie. a script that typically performs the following steps in sequence:The geometry is assembled by gluing brick (rectangular parallelepiped) elements together. The material properties (epsilon and mu) are assumed to be constant within each brick element. Perfect conductors arise as cavities in the domain (no need to glue bricks at conductor locations). The relevant data structure is a dictionary {(i,j,k): [e0,...e11], ...} mapping the brick position, a triplet of integers, (i,j,k) to the 12 edge indices [e0, ...e11]. The sparse matrices A and B are assembled and updated each time a brick element is glued. The input argument "Delta" to the "glue" method is a finite-element, finite difference mixing parameter. Set this "Delta" to 0.9 for improved accuracy (the optimal value may vary between 5/6 and 1). The relevant file is cluster.py with the relevant methods "glue" and "glueFG".
When solving for H, there is no need to modify the sparse matrices A and B, as H satisfies the default natural boundary conditions of the finite element method. When solving for E, however, the amplitude of edges tangent to the conductors must be set to zero (Dirichlet boundary conditions). This can be achieved by setting the diagonal element corresponding to such an edge to arbitrary high value. The module innerfaces.py can be used to extract indices of edges that lie on inner faces. Check also the "assemble" method in curly.py to find out how to modify the A matrix so as to force zero Dirichlet boundary conditions on inner edges.
This is the most cpu time consuming part. Accordingly, Curl3d offers a choice of eigenvalue solvers. For small problems, using SuperLU via the Python interface from Ellipt2d works fine. For larger problems (>10000 variables) consider instead using DistSolve, which interfaces to the Portable, Extension Toolkit for Scientific Computation (PETSc) using an XML-RPC based web service. PETSc offers a large choice of solution methods and preconditioners. Depending on the problem, DistSolve can produce a factor 10 performance improvement over SuperLU. However, DistSolve requires more work to set up. Check method "solve" in curly.py.
Curly3d comes with a simple graphics module ctkvis.py based on the Python Tkinter module. For superior visualization, we recommend the commercial package AVS/Express or any other package reading files in the UCD format. To produce data in UCD format, invoke the method "nodeVector" in ucd.py. By doing so the node values are obtained by interpolating between up to four edges. Alternatively, we also provide an interface to OpenDX (opendx.py).
Methods allowing one to derive the curl, the magnetic field, or the cross product of E and H (to obtain the Poynting flux) are available via edgefilter.py.
Given a driver, run the code simply by invoking
"python curly.py".
where "curly.py" is the driver. You'll need to run your driver in the directory where the Curly3d modules reside, or else define the environment variable PYTHONPATH.
Examples
Send comments to pletzer@pppl.gov
