Comparison of discrete exterior calculus and discrete-dipole approximation for electromagnetic scattering

Journal of Quantitative Spectroscopy & Radiative Transfer ] (]]]]) ]]]–]]]

Contents lists available at ScienceDirect

Journal of Quantitative Spectroscopy & Radiative Transfer journal homepage: www.elsevier.com/locate/jqsrt

Comparison of discrete exterior calculus and discrete-dipole approximation for electromagnetic scattering Jukka Räbinä a,n, Sanna Mönkölä a, Tuomo Rossi a, Antti Penttilä b, Karri Muinonen b,c a b c

Department of Mathematical Information Technology, University of Jyväskylä, P.O. Box 35 (Agora), FI-40014, Finland Department of Physics, University of Helsinki, P.O. Box 64, FI-00014, Finland Finnish Geodetic Institute, Geodeetinrinne 2, P.O. Box 15, FI-02431 Masala, Finland

a r t i c l e i n f o

abstract

Article history: Received 15 October 2013 Received in revised form 31 January 2014 Accepted 10 February 2014

In this study, we perform numerical computations of electromagnetic fields including applications such as light scattering. We present two methods for discretizing the computational domain. One is the discrete-dipole approximation (DDA), which is a well-known technique in the context of light scattering. The other approach is the discrete exterior calculus (DEC) providing the properties and the calculus of differential forms in a natural way at the discretization stage. Numerical experiments show that the DEC provides a promising alternative for solving the general Maxwell equations. & 2014 Elsevier Ltd. All rights reserved.

Keywords: Electromagnetic scattering Discrete exterior calculus Discrete-dipole approximation

1. Introduction We compare two methods for solving the scattering problems numerically. The first method is the discretedipole approximation (DDA), which is a well-known technique in the context of light scattering, see, e.g., [1]. The second approach is based on discrete exterior calculus (DEC), which provides the properties and the calculus of differential forms in a natural way at the discretization stage [2,3]. For electromagnetics the DEC is pioneered by Bossavit and Kettunen [4–6] for antenna applications. Our generalized DEC-formulation for the Maxwell equations [7] works on unstructured grids, and it covers both the classical Yee's FDTD scheme and the Bossavit–Kettunen approach. We concentrate on electromagnetic scattering problems but the methods can be used within other application areas as well. The focus of this paper is

n

Corresponding author. E-mail address: [email protected] (J. Räbinä).

comparing the DEC with the DDA to validate the method and to assess the properties in practice. In particular, we are interested in developing numerical methods tailored for solving large-scale problems efficiently. That is why we pay more attention to the computational efficiency of the methods than to the increase of computing power. The number of computational operations needed for solving the scattering problem with the DDA is of order OðN log NÞ, where N is the number of dipoles used for the space discretization. Assuming that the number of time discretization steps is small compared to the number of space discretization elements N, the computational demand for the DEC method is close to OðNÞ. The rest of this paper is structured as follows. In Section 2, we present the computational model for the light scattering based on the Maxwell equations. We move on with the discretization methods presented in Sections 3 and 4. Section 5 is devoted to numerical experiments and comparison between the methods. Finally, in Section 6, we present concluding remarks.

http://dx.doi.org/10.1016/j.jqsrt.2014.02.011 0022-4073 & 2014 Elsevier Ltd. All rights reserved.

Please cite this article as: Räbinä J, et al. Comparison of discrete exterior calculus and discrete-dipole approximation for electromagnetic scattering. J Quant Spectrosc Radiat Transfer (2014), http://dx.doi.org/10.1016/j.jqsrt.2014.02.011i

J. Räbinä et al. / Journal of Quantitative Spectroscopy & Radiative Transfer ] (]]]]) ]]]–]]]

2

2. Computational model For both discretization methods, the mathematical model of electromagnetic waves is based on the Maxwell equations: ε

∂E ∇  H ¼ sE; ∂t

ð1Þ

μ

∂H þ∇  E ¼ 0; ∂t

ð2Þ

where E and H are the electric and magnetic fields respectively, ε is the electric permittivity, μ is the magnetic permeability, and s is the electric conductivity. The equations presented above describe the time-dependent wave propagation. For considering the time-harmonic wave propagation, we substitute, into Eqs. (1) and (2), the variables E ¼ ReðE expðiωtÞÞ and H ¼ ReðH expðiωtÞÞ, where i is the imaginary unit, such that i2 ¼  1, and ω is the angular frequency. This relationship between timedependent and time-harmonic wave formulations has a crucial role in developing efficient computational methods for time-harmonic waves.

fundamental objects which is called cells (see, e.g., [2]). Basically, a 0-cell is a single point in space, a 1-cell is a differentiable path between two 0-cells, and a p-cell (cell of dimension p) is a connected differentiable p-manifold with a boundary constructed by a finite set of (p 1)-cells. An important property of the differential forms is that a p-form can be integrated on a p-manifold as described in Chapter 7 of Abraham et al. [10]. With this property in mind, the discrete p-form is defined as a set of integration values paired with oriented p-cells of a mesh. For a more detailed description of cells, meshes, and discrete forms we refer the reader to see Hirani's thesis [2]. To discretize the Maxwell equations for the DEC, the vector fields E, B, H, and D must be transferred to the corresponding discrete forms E, B, H, and D. The usual procedure is to define E on the primal edges (1-form), B on the primal faces (2-form), H on the dual edges (1-form), and D on the dual faces (2-form). The values of the discrete forms are computed from vector fields by Z Ei ≔ E  dl; L Zi Bj ≔ B  ds; Sj

3. Discrete exterior calculus and exact controllability approach The discrete exterior calculus provides the properties and the calculus of differential forms (see, e.g., [8]) in a natural way at the discretization stage. For reader who is familiar with the Finite Difference Time-Domain (FDTD) method, we propose to keep that method in mind; in a special case, the DEC is identical to FDTD [5,6]. The advantage of our DEC-based method is the use of arbitrarily shaped well-centered polyhedron elements and alternative Hodge operators. These modifications can significantly increase the accuracy and efficiency of the simulations. Our DEC-based method is developed for solving three-dimensional time-harmonic problems in domains of size 10–100 wavelengths. Next we give a short introduction to the DEC, which is based on the discretization of differential forms. We also present a staggered time discretization scheme and describe an algorithm for utilizing the exact controllability method. A discrete analogue of the smooth differential form is called the discrete differential form (see, e.g., [9,2,3] for details). The spatial discretization is based on a decomposition of the domain into a cell complex, a collection of

Z H j ≔ H  dl; Lnj

Z Di ≔ D  ds; Sni

where Li and Lnj are the primal and dual edge elements, and Sni and Sj are the corresponding dual and primal face elements respectively. It is notable that, with the welldefined mesh and dual mesh, the corresponding dual elements always exist. The positioning of the discrete form values is illustrated in Fig. 1. The space-discretization of the Maxwell equations can be carried out exactly by using the discrete differential forms, and the orientation of a cell complex is an important property in the discretization. The Faraday law is written on the primal elements as n

∂t Bj ¼  ∑ ðd1 Þj;i Ei ; i¼1

where d1 is the incidence matrix defining the neighbouring relations and relative orientations of the primal edges and faces. The incidence matrix d1 corresponds to the curl operator, and it transfers 1-forms to 2-forms. It has nonzero values only, if the current edge Li is included in the

Fig. 1. The E and B fields are defined on the primal elements (left). Similarly the H and D fields are defined on the dual elements (right).

Please cite this article as: Räbinä J, et al. Comparison of discrete exterior calculus and discrete-dipole approximation for electromagnetic scattering. J Quant Spectrosc Radiat Transfer (2014), http://dx.doi.org/10.1016/j.jqsrt.2014.02.011i

J. Räbinä et al. / Journal of Quantitative Spectroscopy & Radiative Transfer ] (]]]]) ]]]–]]]

boundary of the current face Sj. The non-zero entries of the incidence matrix have either the value  1 or the value þ1, depending on the relative orientation of the cell elements [3]. The Ampére law of induction is similarly written on the dual elements as

3

Fig. 2. Staggered time-stepping scheme.

m

J i þ ∂t Di ¼ ∑ ðd1 Þj;i H j ; j¼1

where J is a discrete 2-form on the dual faces. Ji is the free current through Sni . Within the DEC framework, the constitutive relations are replaced by discrete Hodge operators, which are maps from the discrete 1-forms to the discrete 2-forms. Since the discrete forms are presented as vectors, the Hodge operators are logically presented as matrices. The constitutive relations are then D ¼ ⋆ϵE; B ¼ ⋆μH; where ⋆ϵ and ⋆μ are the Hodge operators corresponding to permittivity ϵ and permeability μ respectively. Similarly, the term presenting conductivity can be considered as a Hodge operator ⋆s by J ¼ ⋆sE: The Hodge operators are the only approximated parts in the spatially discretized DEC model: ∂H ¼  ⋆μ  1 d1 E; ∂t   ∂E T ¼ ⋆ϵ  1 d1 H  ⋆sE : ∂t For this reason, there is a lot of research going on concerning the Hodge operators and their accuracy [11,12]. The construction of the mesh and the corresponding dual mesh has an important role on accuracy of the Hodge operator. In principle, polyhedra of different shapes can be used to construct a mesh for the DEC method. Typically, the elements are wanted to be well-centered to guarantee orthogonal duality. The orthogonality of the primal and dual elements implies diagonal Hodge operators providing a significant saving in computing time. This can be obtained by using the weighted circumcenters for constructing the dual mesh. Physical applications are defined in unbounded, that is, exterior, domains. To solve the problems numerically, the domain is truncated to a finite computational domain by boundary conditions. In practice, an artificial boundary condition, that is an approximation of the Sommerfeldtype radiation condition, is set on the exterior artificial boundary, and it ensures that the solution approximates the restriction of the solution in the original unbounded region. With the DEC, we use the absorbing layers based on Berenger's original idea of a perfectly matched layer (PML) [13]. Since we are considering the time-dependent wave equation, we must proceed with the simulation of time evolution. The time discretization is established with a staggered leapfrog-type iteration in which the values for the variables E and H are computed at separate time instances as shown in Fig. 2. The time discretization scheme is second-order accurate. With a diagonal Hodge

operator the method is explicit, i.e., the values for each time step are determined from the values of the previous time step. The computations can be performed in an efficient manner, since only matrix–vector products are needed in time-dependent simulation. The stability of the leapfrog manner iteration depends on the length of the time step Δt. If the time step is too long, it will cause E and H components bouncing between the negative and positive values and the absolute values will increase without limit as the simulation progresses. That is why the scheme needs to satisfy the Courant– Friedrichs–Lewy (CFL) condition [14], which guarantees the numerical stability but limits the length of the time step. Instead of setting a constant time step size Δt, it is often computationally more efficient to assign each element by its own, optimized time step, as done by Stern et al. [3]. This asynchronous time stepping scheme is based on the method of asynchronous variational integrators (AVI), which was earlier applied for problems in elastodynamics by Lew et al. [15]. Our time discretization scheme differs slightly from the AVI. We use synchronized non-uniform time stepping, where each time step size is expressed by Δt=3n , n A N. Thus, we get a simple non-uniform time stepping scheme, which conserves some good convergence properties of staggered leapfrog discretization. Within the framework defined above, both timeharmonic and time-dependent electromagnetic problems can be considered. With this model we can, for instance, solve a time-harmonic scattering problem and then simulate how the solution evolves in time. Since our focus is on electromagnetic scattering problems having sinusoidal time variation, our aim is to get time-harmonic solutions. Essentially, the time-harmonic solution can be reached by simple time integration. Nevertheless, such an asymptotic approach suffers from poor convergence in particular for large wave numbers and complicated domains. Thus, we accelerate the convergence rate by using the exact controllability technique, first introduced within scalar-valued wave equations in [16]. Recently, it has been applied to acoustics and elastodynamics (see, e.g., [17,18]). The theoretical background of combining the exact controllability method with discrete exterior calculus for the generalized time-periodic wave equations is presented by Pauly and Rossi [7]. The key point of the exact controllability technique is searching for the periodic solution of the time-dependent wave equation as a leastsquares optimization problem. In principle, the timeharmonic solution can be found by minimizing the difference between the initial conditions of the time-dependent problem and the corresponding variables after one time period. A natural quadratic error functional is the squared energy norm of the system measuring how close the solution is to the corresponding time-harmonic solution. The energy norm is a weighted L2-norm, implying that the

Please cite this article as: Räbinä J, et al. Comparison of discrete exterior calculus and discrete-dipole approximation for electromagnetic scattering. J Quant Spectrosc Radiat Transfer (2014), http://dx.doi.org/10.1016/j.jqsrt.2014.02.011i

J. Räbinä et al. / Journal of Quantitative Spectroscopy & Radiative Transfer ] (]]]]) ]]]–]]]

4

discrete quadratic functional we minimize is spanned by a diagonal mass matrix. In practice, the minimization is performed by the conjugate gradient (CG) method operating in a purely L2-type Hilbert space. Thus, the minimization problem is well-conditioned, and the number of CG iterations does not grow with the mesh refinement.

4. Discrete-dipole approximation The discrete-dipole approximation (DDA) is a wellestablished technique and its details and formulation are already covered in many publications. We recommend the interested reader to refer to, i.e., Yurkin and Hoekstra [1] for the general description of the method. In short, the DDA is a discretized volume-integration method to solve the integral equation of the electric field inside the scatterer, i.e., in domain V, Z 3 EðrÞ ¼ Einc ðrÞ þ d r 0 Gðr; r0 Þχðr 0 Þ Eðr 0 Þ; ð3Þ

The possible errors in the DDA method can arise from two sources. First, the target geometry is discretized and this can introduce geometrical errors. Second, the polarizability description might not be exact, especially for the dipoles at the boundary of the target. Both these error sources can be made practically negligible by decreasing the dipole size. It seems that the theoretical results in this field are difficult [19], but there are many numerical studies (e.g., [20–22]) that show good convergence in practice with decreasing dipole size when DDA is applied to targets where jr  1j o 2 with r as the complex refractive index of the target. There are at least two excellent open-source implementations of the DDA method with tested accuracy [23], and we choose here to use the ADDA code, version 1.2 [24]. Several iterative methods for solving the DDA equations are implemented in ADDA, and we employ here the biconjugate gradient stabilized method and the quasi-minimal residual method. For the dipole polarizability, the lattice dispersion relation is used.

V

inc

where E ðrÞ and EðrÞ are respectively the incident and the total electric field at location r, χ is the susceptibility of the medium, and G is the free space dyadic Green's function [1]. If the domain V is divided into subvolumes Vi, and the field and the susceptibility are considered constant within each Vi, the discretized form of Eq. (3) can be written as 1 Einc i ¼ α i Pi  ∑ G ij Pj ; iaj

ð4Þ

¼ V i χ i Ei , α i are the where dipole polarizations Pi ¼ α i Eexc i polarizability tensors, and Eexc is the exciting electric fields. i Different formulations, such as Clausius–Mossotti, can be used to calculate the polarizabilities in each subvolume, i. e., dipole. The form of Eq. (4) will lead to a group of linear equations that are solved in the DDA method. The linear system in DDA is usually quite large since it has 3N unknowns where N is the number of dipoles. Direct solution becomes impossible, and iterative methods are used. Although it is not required in the DDA that the dipoles should lie in a rectangular grid, in practice this is often desired. With rectangular grid the discrete fast Fourier transform can be used to speed up the matrix– vector multiplications in the iterative solver.

5. Numerical experiments In this section we compare the efficiency of the DEC and the DDA by performing electromagnetic scattering simulations. By efficiency of a method we refer to the ratio between the accuracy and the computational work required to achieve that level of accuracy. The accuracy of the discrete solution depends on space discretization, stopping criterion of the iterative solving procedure and approximation of the geometry. In addition, with the DECbased approach there are additional sources of inaccuracy. These are the approximation of the radiation condition, constructed as an absorbing layer, and the time integration performed by the leapfrog scheme. We compare the methods by solving electromagnetic scattering problems, given by the model (1) and (2), with ε ¼ 2:5575, μ ¼ 1:0, and s ¼ 0:16 implying the refractive index m ¼ 1:6 þ 0:05i. The incident wave is chosen to be a fully polarized plane wave of wavelength λ ¼ 2π propagating in the direction of the positive x-axis. Analytical solutions for the scattering problems of this kind exist only for some simple geometries, like a sphere. That is why we have chosen to use a sphere of diameter 5λ

Fig. 3. The DEC discretizations by polygons of different shapes (tetrahedron, prism, cube, octahedron). For clarity, these figures are made for 4 elements per wavelength, but 10–30 elements per wavelength are used in the simulations. An absorbing layer of thickness 1:5λ is surrounding the scattering objects. (For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this paper.)

Please cite this article as: Räbinä J, et al. Comparison of discrete exterior calculus and discrete-dipole approximation for electromagnetic scattering. J Quant Spectrosc Radiat Transfer (2014), http://dx.doi.org/10.1016/j.jqsrt.2014.02.011i

J. Räbinä et al. / Journal of Quantitative Spectroscopy & Radiative Transfer ] (]]]]) ]]]–]]]

as an obstacle in the first numerical experiment. In this case, the accuracy of the methods is validated by comparing the solution with the known analytical solution available by the Mie theory. The other two obstacles, which we use, are a cube and a torus. Both the edge length of the cube and the outer diameter of the torus are set to be 5λ. The inner diameter of the torus is set to be 3λ, which means that the thickness of the torus ring is 2λ. The incident wave propagates orthogonal to one of the cubic faces and through the torus hole. All the obstacles are illustrated by yellow colour in Fig. 3. In addition, Fig. 3 illustrates the computational domain and the space discretization, which are used with DECbased method. The absorbing layer of thickness 1:5λ, based on Berenger's original idea of the perfectly matched layer [13], is surrounding the obstacles. Berenger's PML is designed to absorb the outgoing wave and leave the orthogonal wave unaffected. Thus, the idea of separating the outgoing wave components from the orthogonal wave components is imported in the DEC system. The separation of components requires specific mesh design where primal and dual edges are either parallel or orthogonal to the absorption direction. The absorption of the outgoing wave is accomplished by changing electric and magnetic conductivities sx and snx by the depth in layer x. Starting from the zero at the inner level x¼0, the conductivities are linearly increased by the relation sx sn ¼ x ¼ βx; ϵ0 μ0 In these simulations, constant β ¼ 0:7 is used. The absorption directions of each problem are illustrated in Fig. 4. The volumes with several absorption directions are treated in a special manner. The meshes are constructed by combining parts of structured non-regular grids, where elements are wellcentered. Between the structured parts there are unstructured elements which are then optimized by the Hodge optimization [11]. In this optimization, the element positions and weights are corrected to have more stable Hodge properties. The construction from structured parts is relatively fast compared to fully unstructured mesh construction. Also, element sizes and properties can be easily controlled.

5

Although the resolution of the space discretization presented in Fig. 3 is 4 elements per wavelength, in the computations we use, in principle, 10, 15, 20, and 30 elements per wavelength. In practice, there is some variation in mesh step sizes in the unstructured parts of the mesh. Element sizes are adjusted to fit the local wavelength. The elements of the DDA system are dipoles placed on a regular grid. In that case, the computational domain is the minimal rectangular parallelepiped that encompasses the scatterer. No absorbing layers are needed like in the DEC. The resolutions of the space discretization used with the DDA are 10, 20, and 30 dipoles per wavelength. The cube and the sphere are discretized by the systems of 80  80  80, 160  160  160 and 240  240  240 dipoles, whereas the torus is discretized by the systems of 16  80  80, 32  160 160, and 48  240  240 dipoles respectively. For the simulations with obstacles of non-spherical shapes, the reference solution is computed by the DDA method with 80 dipoles per wavelength. In this case, the number of dipoles is 640  640  640 for the cube and 128  640  640 for the torus. This dense space discretization gives a highly accurate solution. In the case of the cube, the phase function error is known to be less than 1 percent at every measurement angle. The phase functions of each simulations are illustrated in Fig. 5. In practice, the accuracy of both methods is considered as the relative error of the phase function, i.e., by comparing the simulated phase functions Ps to the reference phase function P0 by the formula Rπ jP s ðθÞ P 0 ðθÞj sin θ dθ δP s ¼ 0 R π : 0 P 0 ðθÞ sin θ dθ The far-field scattering data is obtained from the DEC result by applying a powerful and flexible near-field to farfield transformation of Taflove and Umashankar [25]. The near-field data was collected from a surface outside the scatterer as illustrated in Fig. 4. The transformation from the discrete differential forms to the field vectors is always a source of error in the DEC-based method. To minimize this error, we have constructed the mesh for data collection such that the distance of 0:2λ is used between the near-field surface and the scatterer. Thus, the full dimensions of the DEC domain are up to 8:4λ.

Fig. 4. The scattering fields of the DEC simulation are presented by background colours, where red, green, and blue components represent the x-, y- and zcomponents of the electric field respectively. The white line represents the surface, where the near-field data is collected. The black lines illustrate the absorption directions of the PML layers. The cross sections in the xy-plane are illustrated. (For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this paper.)

Please cite this article as: Räbinä J, et al. Comparison of discrete exterior calculus and discrete-dipole approximation for electromagnetic scattering. J Quant Spectrosc Radiat Transfer (2014), http://dx.doi.org/10.1016/j.jqsrt.2014.02.011i

6

J. Räbinä et al. / Journal of Quantitative Spectroscopy & Radiative Transfer ] (]]]]) ]]]–]]]

Fig. 5. The scattering phase functions of the spherical, cubic and toroidal objects are computed with 30 discretization elements per wavelength. The reference phase function of the spherical object is obtained by the exact Mie theory. The highly accurate DDA solutions with 80 elements per wavelength are used as the reference solutions for the simulation results with the cube and the torus.

Fig. 6. Relative error of the phase function with respect to the simulation time in seconds.

The residual of the DEC-based algorithm defines, at each iteration, how far the solution is from a periodic one. The CG iterations are continued until the relative norm of the residual is smaller than 10  4. In these tests, only 30– 35 iterations were needed. The DEC simulations were run by the algorithm developed at the University of Jyväskylä. The algorithm is implemented in the Cþ þ programming language. The DDA simulations were run by the ADDA DDA v. 1.2 code. The stopping criterion for the iterative solver in ADDA, concerning the relative norm of the residual, was set to 10  8. With the cube shape, the Bi-CG stabilized method was used, but with the sphere and the torus the iteration was successful using the default method

of quasi-minimal residual. All the computations have been carried out on an Intel Xeon E5-2670 processor at 2.6 GHz. As the efficiency of a numerical simulation method corresponds to the level of accuracy in a given amount of simulation time, the relative error of the phase function is presented as a function of time in Fig. 6. It can be seen that the computational cost is of the same order of magnitude for both methods. Further, it is possible to improve the accuracy at the expense of computational time with both methods by using denser space discretization. As we take a closer look, we can see that there are small differences in the performance of the two methods. The accuracy is slightly better for a cubic obstacle with the DDA than with the DEC. That is because of the nature of the discretization with the DDA. The cube fills optimally the computational domain, and the full dipole grid is used efficiently. There is no need to use any additional dipoles outside the obstacle, which gives savings in computing time compared to the DEC. It is worth mentioning that the DEC discretization in these experiments is constructed using a tetrahedral honeycomb tessellation. Another option would be using cubic elements like in the classical FDTD method. On the other hand, the reference solutions are solved by the DDA. Thus, it is natural that the DDA solution converges towards the more accurate DDA solution. For both methods, it is challenging to approximate accurately the geometries of curvilinear shapes, such as the boundary of a sphere or a torus. In particular with the

Please cite this article as: Räbinä J, et al. Comparison of discrete exterior calculus and discrete-dipole approximation for electromagnetic scattering. J Quant Spectrosc Radiat Transfer (2014), http://dx.doi.org/10.1016/j.jqsrt.2014.02.011i

J. Räbinä et al. / Journal of Quantitative Spectroscopy & Radiative Transfer ] (]]]]) ]]]–]]]

DDA, the structured discretization does not fit well on modelling curved boundaries. In the DEC, there are two ways to treat the curved boundaries. The first method is to build a mesh where the curved boundaries are approximated by cells (as we have carried out here, see Fig. 3). The second way is to share the cells between the objects with different material parameters. Sharing can be established easily by setting a cell material parameter as the weighted average of the object material parameters. We can see that both the DDA and the DEC methods solve the scattering problem by a sphere with about the same efficiency. Further, we can observe that the accuracy is somewhat better with the DEC for a torus. The toroidal object was incorporated into the present study due to its sparse spatial characteristics. In the DEC discretization, the unstructured mesh with changing mesh step size is used to better approximate the complex geometry. This is important in the case of sparse objects. On the other hand, the toroidal object is the smallest of the objects under present study. Hence, in this particular test case, a considerably large absorbing layer is needed for the DEC. Therefore, a relatively large amount of computing time is used for processing with the absorbing layer. Consequently, which of the two methods is more efficient depends on the shape and the size of the obstacle. The DEC seems to be a promising method in the light of accuracy and computational cost especially for nonconvex obstacles. 6. Conclusions We compared two methods, the discrete exterior calculus (DEC) and the discrete-dipole approximation (DDA), for solving electromagnetic scattering problems for three different obstacles. The computational efficiencies are of the same order of magnitude with both methods. We can also conclude that both of the approaches have their own strengths and weaknesses, depending on the geometry of the scatterer and the details of the discretization. For instance, curved shapes, such as a torus, can be approximated with high accuracy within the DEC framework. On the other hand, the DDA works better with cubic obstacles. Combining the exact controllability method with the DEC provides a promising alternative for solving the general Maxwell equations. Combining cubical elements with arbitrary shaped well-centered elements provides geometric flexibility, still maintaining the computational efficiency. Also the implemented asynchronous timestepping makes the simulation efficient at every point of the geometry. Benefits of the DEC-based approach include also implementability for parallel computers and multicore processors, testing of which is left as a future work. A lot of improvement can still be done also by developing the mesh generation that serves the accuracy of the Hodge operator. Acknowledgements The computations presented have been made using the computing resources provided by the CSC — IT center for

7

science, owned by the Finnish Ministry of Education and Culture. This project has been funded by the Academy of Finland, Grants 259925, 257966, and 260076. References [1] Yurkin MA, Hoekstra AG. The discrete dipole approximation: an overview and recent developments. J Quant Spectrosc Radiat Transf 2007;106(1–3):558–89. [2] Hirani AN. Discrete exterior calculus [Ph.D. thesis]. California Institute of Technology; 2003. [3] Stern A, Tong Y, Desbrun M, Marsden JE. Geometric computational electrodynamics with variational integrators and discrete differential forms, preprint, arXiv:0707.4470v3 [math.NA], 2009. [4] Bossavit A. In: Computational electromagnetismvariational formulations, complementarity, edge elements. San Diego: Academic Press; 1998. [5] Bossavit A, Kettunen L. Yee-like schemes on a tetrahedral mesh, with diagonal lumping. Int J Numer Model 1999;12(1–2):129–42. [6] Bossavit A, Kettunen L. Yee-like schemes on staggered cellular grids: a synthesis between fit and FEM approaches. IEEE Trans Magn 2000;36(4):861–7. [7] Pauly D, Rossi T. Theoretical considerations on the computation of generalized time-periodic waves. Adv Math Sci Appl 2011;21(1): 105–31. [8] Schutz BF. Geometrical methods of mathematical physics. Cambridge, UK: Cambridge University Press; 1980. [9] Desbrun M, Hirani AN, Leok M, Marsden JE. Discrete exterior calculus, preprint, arXiv:math/0508341v2 [math.DG], 2005. [10] Abraham R, Marsden JE, Ratiu T. Manifolds, tensor analysis, and applications (applied mathematical sciences). New York: SpringerVerlag; 1988. [11] Mullen P, Memari P, de Goes F, Desbrun M. HOT: hodge-optimized triangulations. ACM Trans Graph 2011;30(4):103. [12] Hirani AN, Kalyanaraman K, VanderZee EB. Delaunay Hodge star, preprint, arXiv:0707.4470v3 [math.NA] (2012). [13] Berenger J-P. A perfectly matched layer for the absorption of electromagnetic waves. J Comput Phys 1994;114(2):185–200. [14] Courant R, Friedrichs K, Lewy H. On the partial difference equations of mathematical physics. IBM J 1967:215–34 (English translation of the 1928 German original). [15] Lew A, Marsden JE, Ortiz M, West M. Asynchronous variational integrators. Arch Rational Mech Anal 2003;167(2):85–146. [16] Bristeau MO, Glowinski R, Périaux J. Controllability methods for the computation of time-periodic solutions; application to scattering. J Comput Phys 1998;147(2):265–92. [17] Kähkönen S, Glowinski R, Rossi T, Mäkinen R. Solution of timeperiodic wave equation using mixed finite-elements and controllability techniques. J Comput Acoust 2011;19(4):335–52. [18] Mönkölä S, Heikkola E, Pennanen A, Rossi T. Time-harmonic elasticity with controllability and higher order discretization methods. J Comput Phys 2008;227(11):5513–34. [19] Yurkin MA, Maltsev VP, Hoekstra AG. Convergence of the discrete dipole approximation. I. Theoretical analysis. J Opt Soc Am A 2006;23(10):2578–91. http://dx.doi.org/10.1364/JOSAA.23.002578. [20] Draine B, Flatau P. The discrete dipole approximation for scattering calculations. J Opt Soc Am A 1994;11:1491–9. http://dx.doi.org/ 10.1364/JOSAA.11.001491. [21] Yurkin MA, Maltsev VP, Hoekstra AG. Convergence of the discrete dipole approximation. II. An extrapolation technique to increase the accuracy. J Opt Soc Am A 2006;23(10):2592–601. http://dx.doi.org/ 10.1364/JOSAA.23.002592. [22] Yurkin MA, Kahnert M. Light scattering by a cube: accuracy limits of the discrete dipole approximation and the T-matrix method. J Quant Spectrosc Radiat Transf 2013;123:176–83 http://dx.doi.org/10.1016/ j.jqsrt.2012.10.001. [23] Penttilä A, Zubko E, Lumme K, Muinonen K, Yurkin M, Draine B, et al. Comparison between discrete dipole and exact techniques. J Quant Spectrosc Radiat Transf 2007;106:417–36. http://dx.doi.org/10.1016/ j.jqsrt.2007.01.026. [24] Yurkin M, Hoekstra A. The discrete-dipole-approximation code ADDA: capabilities and known limitations. J Quant Spectrosc Radiat Transf 2011;112(13):2234–47. http://dx.doi.org/ 10.1016/j.jqsrt.2011.01.031. [25] Taflove A, Umashankar K. Radar cross section of general threedimensional scatterers. IEEE Trans Electromagn Compat 1983;25(4): 433–40.

Please cite this article as: Räbinä J, et al. Comparison of discrete exterior calculus and discrete-dipole approximation for electromagnetic scattering. J Quant Spectrosc Radiat Transfer (2014), http://dx.doi.org/10.1016/j.jqsrt.2014.02.011i