We present the discontinuous Galerkin method combined with a low-storage Runge-Kutta method as an accurate and efficient way to numerically solve the time-dependent Maxwell’s equations. We investigate the numerical scheme in the context of modeling scattered electromagnetic wave’s propagation through human eye’s structures

Araújo, Adérito

Barbeiro, Sílvia

Bernardes, Rui

Caramelo, Francisco

Correia, António L.

Morgado, Miguel

Pinto, Luís

Serranho, Pedro

Silva, Ana S.C.

Estudo Geral

Proceedings of the 13th International Conference
on Computational and Mathematical Methods
in Science and Engineering, CMMSE 2013
24–27 June, 2013.
Numerical solution of time-dependent Maxwell’s equations
for modeling scattered electromagnetic wave’s propagation
Ade´rito Arau´jo1, S´ılvia Barbeiro1, Lu´ıs Pinto1, Francisco Caramelo2,
Anto´nio L. Correia2, Miguel Morgado2,5, Pedro Serranho2,4, Ana S´ılvia F.
C. Silva5 and Rui Bernardes2,3
1 CMUC, Department of Mathematics, University of Coimbra, Portugal
2 IBILI, Faculty of Medicine, University of Coimbra, Portugal
3 ABILI, Coimbra, Portugal
4 Mathematics Section, Department of Science and Technology, Open University, Portugal
5 Department of Physics, Faculty of Science and Technology, University of Coimbra,
Portugal
emails: alma@mat.uc.pt, silvia@mat.uc.pt, luisp@mat.uc.pt,
fcaramelo@fmed.uc.pt, acorreia@aibili.pt, mmorgado@ibili.uc.pt,
Pedro.Serranho@uab.pt, anascsilva@gmail.com, rmbernardes@fmed.uc.pt
Abstract
We present the discontinuous Galerkin method combined with a low-storage Runge-
Kutta method as an accurate and e cient way to numerically solve the time-dependent
Maxwell’s equations. We investigate the numerical scheme in the context of modeling
scattered electromagnetic wave’s propagation through human eye’s structures.
Key words: Maxwell’s equations, discontinuous Galerkin method, low-storage Runge-
Kutta method, anisotropic permittivity tensor
1 Introduction
The human retina is a complex structure in the eye that is responsible for the sense of
vision. It is part of the central nervous system, constituted by layers of neurons intercon-
nected through synapses [7]. There are ten layers in total; one of them is constituted by
c CMMSE ISBN: 978-84-616-2723-3Page 121 of 1797
Numerical solution of time-dependent Maxwell’s equations
photosensitive neurons. There are a number of eye-related pathologies that can be identified
by analysis of these retinal layers in detail [7]. All these pathologies can be diagnosed more
conclusively with the help of the increasingly popular optical imaging technique – optical
coherence tomography (OCT) [3], [15]. In fact, previous studies have established a link
between changes in the blood-retina barrier and in optical properties of the retina [1], [2]
which can be identified by this exam.
Physically, OCT it is based in low coherence interferometry. This technique uses an
electromagnetic wave with a low coherence length. In order to better understand the in-
formation carried in an optical coherence tomography, it is crucial to study in detail the
behaviour of the electromagnetic wave as it travels through the sample. Several di↵erent
models have been developed to describe the interactions of the electromagnetic field with
biological structures. The first models were based on single-scattering theory [14], which
is restricted to superficial layers of highly scattering tissue in which only single scattering
occurs. Simulating the full complexity of the retina, in particular the variation of the size
and shape of each structure, distance between them and the respective refractive indexes,
requires a more rigorous approach that can be achieved by solving Maxwell’s equations.
In this work we discuss the numerical discretization of the time-dependent Maxwell’s
equations. We use the discontinuous Galerkin (DG) method for the integration in space and
a low-storage Runge-Kutta method for the integration in time. In the model we consider
anisotropic permittivity tensors which arise naturally in our application of interest. We
illustrate the performance of the method with some numerical experiments.
2 Maxwell’s Equations
We shall consider the time domain Maxwell’s equations in the two-dimensional transverse
electric (TE) mode. For this case, and assuming no conductivity e↵ects, the equations in
the non-dimensional form are
✏
@Ex
@t
=
@Hz
@y
(1)
✏
@Ey
@t
=  @H
z
@x
(2)
µ
@Hz
@t
=  @E
y
@x
+
@Ex
@y
, in ⌦⇥ (0, T ], (3)
where E = (Ex, Ey) and Hz represent the electric and magnetic fields, respectively, ✏
represents the relative permittivity of the medium and µ is the permeability of the medium
and ⌦ is a two-dimensional domain. The nondimensionalized variables in (1)-(3) are related
to the physical variables in the following way
x˜
L
= x,
y˜
L
= y,
ct˜
L
= t, E = Z 10 E, H
z = Hz,
c CMMSE ISBN: 978-84-616-2723-3Page 122 of 1797
A. Arau´jo, S. Barbeiro, L. Pinto, R. Bernardes et al.
where L is a reference length, c is the speed of light in free space, and Z0 =
p
µ0/✏0 is
the free-space impedance. The set of equations (1)-(3) must be complemented by proper
boundary conditions, for instance, the perfect electric boundary condition (PEC) [5]
⌘ ⇥ E = 0, on @⌦, (4)
or the Silver-Mu¨ller absorbing boundary condition [9]
⌘ ⇥ E =
r
µ
✏
⌘ ⇥ (Hz ⇥ ⌘), on @⌦. (5)
In both cases ⌘ denotes the unit outward normal. To complete the model, initial conditions
E0 = E(0) and H0 = H(0), in ⌦,
must also be considered.
3 The scattered-field formulation
For the investigation of scattering problems in linear materials, we exploit the linearity of
the Maxwell’s equations (1)-(3) and separate the fields (E, H) into incident (Ei, H i) and
scattered components (Es, Hs), i.e.,
E = Es + Ei and H = Hs +H i. (6)
Assuming that the incident field is also a solution of the Maxwell’s equations we obtain,
after some manipulation the scattered field formulation as in [16],
✏
@Ex,s
@t
=
@Hz,s
@y
+ P (7)
✏
@Ey,s
@t
=  @H
z,s
@x
+Q (8)
µ
@Hz,s
@t
=  @E
y,s
@x
+
@Ex,s
@y
+R, in ⌦⇥ (0, T ] (9)
with the source terms
P (x, y, t) = (✏i   ✏)@E
x,i
@t
,
Q(x, y, t) = (✏i   ✏)@E
y,i
@t
,
R(x, y, t) = (µi   µ)@H
z,i
@t
,
where ✏i and µi represent the relative permittivity and permeability of the medium in which
the incident field propagates. In our research, for simplicity, we adopt the pure scattered
field formalism in opposition to the total field/scattered field formalism [16].
c CMMSE ISBN: 978-84-616-2723-3Page 123 of 1797
Numerical solution of time-dependent Maxwell’s equations
4 Numerical method
The DG method was first introduced in [11]. In particular, the nodal formulation described
in [6] has gained notorious popularity in recent years and it has been extensively used in
electromagnetic problems since the first application of the method to Maxwell’s equations
in 2002 (see [5]). This is the method that we have chosen to use. In opposition to the
traditional finite di↵erence (FD) time domain methods [16] based on the Yee’s scheme [17]
the DG method is a high-order accurate method that can easily handle complex geometries.
Moreover, local refinement strategies can easily be incorporate due to the ability of the
method to deal with irregular meshes with hanging nodes and local spaces of di↵erent orders.
When compared to finite element methods [12], the DG method presents the advantages of
avoiding the solution of linear systems when explicit time integrators are employed and the
suitability for parallel implementation on modern multi-graphics processing units (GPUs)
[6].
As starting point of our work we use the MatLab codes for the DG method [6]. The
software includes many attributes that we intent to exploit in future work such as three-
dimensional routines and nonconforming triangulations. For now, we are concerned with
two-dimensional problems and we have restricted our attention to conform triangular ele-
ments. Our main focus has been on the implementation of additional features needed for
the kind of problems which are in our objectives. For instance, the retinal nerve fiber layer
of the retina is a birefringent medium, meaning that tensorial permittivity ✏ needs to be
taken into account. We addressed this issue using the numerical scheme presented in [8].
The procedure, based on the so called upwind flux, is fairly simple and can be easily in-
corporated in the algorithm. Nevertheless, it seems to be e cient and robust for a wide
range of problems, since the two by two matrix that represents the tensor ✏ only needs
to be invertible. In [8] the authors claim that the method, which is an extension of the
basic two-dimensional formulation for isotropic materials to allow anisotropic permittivity
tensors, retains the convergence characteristics of the original method. The validation tests
that we present in the next section corroborate those findings.
Another aspect of our study concerns the integration in time. The time integration
scheme used in the original algorithm [6] is a fourth order, five stage low-storage Runge-
Kutta (LSRK) method. Such type of schemes are very popular in this context as they retain
the qualities of the original Runge-Kutta schemes while decreasing the memory consumption
significatively. However, we recall that explicit methods like the LSRK scheme are subject
to the CFL stability condition. In practice, this means that the time step size is proportional
to the smallest elements of the spatial mesh [6]. For some problems this condition can have
a strong impact on the e ciency of the algorithm. In the future we intent to deal with this
issue using local time-stepping strategies [6] or locally implicit time-schemes [4]. For now, we
mitigate this restriction implementing the improved fourth order, 14-stage LSRK presented
in [10]. With this scheme the authors report a speed improvement of about 40%  50% in
c CMMSE ISBN: 978-84-616-2723-3Page 124 of 1797
A. Arau´jo, S. Barbeiro, L. Pinto, R. Bernardes et al.
relation to the previous fourth order, five stage LSRK method. This gain is consequence of
an improved stability region possible by a suitable choice of the coe cients of the scheme.
This is a very attractive approach, with very few changes in the code and no additional cost
in memory or accuracy. The validation tests that we present in the next section indicate
that these conclusions are also valid for the tensorial case.
An important aspect in computational electromagnetic problems is the implementa-
tion of absorbing boundary conditions. In our work we have implemented the Silver-
Mu¨ller boundary condition (5) and the well established and more e↵ective uniaxial perfectly
matched layer (UPML) [13]. The UPML can be incorporated in the DG method without
any major modification and so far has shown to be an e cient approach.
5 Numerical results
We consider initial conditions, boundary conditions and functions P (t), Q(t) and R(t), such
that the system of equations (7)-(9) has the solution
Ex = t2(  cos(⇡x)  1)(1
3
y3   y) (10)
Ey = 0 (11)
Hz =
1
3
t3(  cos(⇡x)  1)(y2   1) (12)
with ✏ = 4x2c + y
2
c + 1, for the results in Tables 1 and 2, and
✏ =

4x2c + y
2
c + 1 |xc + yc|
|xc + yc| x2c + 2
 
. (13)
for the results in Tables 3-6. In both cases we set µ = 1, ⌦ = [ 1, 1]2 and T = 5. By xc
and yc we represent the centroid of the triangles. Note that in the implementation ✏ must
be constant in each triangle.
hmax kEx   ExhkL2 Rate kEy   EyhkL2 Rate kHz  HzhkL2 Rate
1.4142e-01 4.1644e-01 2.0580 4.1702e-01 2.0068 3.5908e-01 2.0297
7.0711e-02 2.4017e-02 2.0265 2.5820e-02 2.0002 2.1537e-02 2.0101
3.5355e-02 1.4469e-03 - 1.6133e-03 - 1.3273e-03 -
Table 1: Convergence rate in the L2 norm for polynomial approximation of order 1.
From the analysis of Tables 1-4 we observe that the tensorial scheme preserves the
optimal rate of convergence of the upwind flux in the scalar case O(hN+1), where N denotes
the order of the polynomial involved in the approximation. Note also that the order of the
c CMMSE ISBN: 978-84-616-2723-3Page 125 of 1797
Numerical solution of time-dependent Maxwell’s equations
hmax kEx   ExhkL2 Rate kEy   EyhkL2 Rate kHz  HzhkL2 Rate
1.4142e-01 8.6363e-04 3.0268 1.0196e-03 3.0124 1.7076e-03 2.9941
7.0711e-02 1.3002e-05 3.0120 1.5660e-05 3.0070 2.6901e-05 2.9978
3.5355e-02 1.9981e-07 - 2.4231e-07 - 4.2159e-07 -
Table 2: Convergence rate in the L2 norm for polynomial approximation of order 2.
hmax kEx   ExhkL2 Rate kEy   EyhkL2 Rate kHz  HzhkL2 Rate
1.4142e-01 5.4515e-01 2.0461 6.2046e-01 1.9949 3.9240e-01 2.0202
7.0711e-02 3.1963e-02 2.0197 3.9055e-02 1.9964 2.3847e-02 2.0067
3.5355e-02 1.9438e-03 - 2.4532e-03 - 1.4767e-03 -
Table 3: Convergence rate in the L2 norm for polynomial approximation of order 1.
hmax kEx   ExhkL2 Rate kEy   EyhkL2 Rate kHz  HzhkL2 Rate
1.4142e-01 8.1456e-04 2.9993 9.4071e-04 3.0184 1.7020e-03 2.9940
7.0711e-02 1.2739e-05 2.9996 1.4328e-05 3.0091 2.6817e-05 2.9979
3.5355e-02 1.9916e-07 - 2.2106e-07 - 4.2022e-07 -
Table 4: Convergence rate in the L2 norm for polynomial approximation of order 2.
hmax kEx   EhxkL2 Time
original fourth order, 5-stage LSRK 7.0711e-02 3.1963e-02 69.25s
improved fourth order, 14-stage LSRK 7.0711e-02 3.1963e-02 38.35s (44%)
original fourth order, 5-stage LSRK 3.5355e-02 1.9438e-03 514.26s
improved fourth order, 14-stage LSRK 3.5355e-02 1.9438e-03 295.99s (42%)
Table 5: Comparison of the two LSTR methods for polynomial approximation of order 1.
hmax kEx   EhxkL2 Time
original fourth order, 5-stage LSRK 7.0711e-02 1.2739e-05 132.43s
improved fourth order, 14-stage LSRK 7.0711e-02 1.2739e-05 79.71s (39%)
original fourth order, 5-stage LSRK 3.5355e-02 1.9916e-07 1259.94s
improved fourth order, 14-stage LSRK 3.5355e-02 1.9916e-07 745.40s (40%)
Table 6: Comparison of the two LSTR methods for polynomial approximation of order 2.
errors is about the same when we compare the problem with tensorial permittivity and the
problem with scalar permittivity.
In Tables 5 and 6 we compare the two di↵erent implementations for the time integration.
The error is exactly the same for both methods and the running time decreases about 40%.
c CMMSE ISBN: 978-84-616-2723-3Page 126 of 1797
A. Arau´jo, S. Barbeiro, L. Pinto, R. Bernardes et al.
Even though we only show the results for Ex, in our experiments we found that the error
of the other field components, Ey and Hz, present the same type of behaviour.
To illustrate the advantages of the DG method over the FD method we examine the
scattering by a dielectric cylinder with a relative permittivity of ✏ = 5. This scatterer with
radius r = 0.5 and centered at the origin is excited by a Gaussian plane wave propagating
along the y direction. We measured the magnetic filed Hz at the point (x, y) = ( 1, 0).
In Figure 1 we compare the results obtained with both methods. The DG solution was
obtained using third-order polynomials and a local refined mesh formed by 920 triangles.
This means that the number of unknowns is 9.2⇥ 103. The FD solution was obtained using
a uniform mesh with h = 1.25⇥ 10 2 implying 4.0⇥ 104 unknowns. Our experiments show
that the DG method is more accurate and e cient. The error, in the Euclidean norm, is
1.5827⇥ 10 1 for the DG method and 3.6668⇥ 10 1 for the FD method.
0 1 2 3 4 5 6 7 8 9
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
Time
Figure 1: Error curves for the DG solution (dash line) and FD solution (solid line).
Acknowledgements
This work has been partially supported by Centro de Matema´tica da Universidade de
Coimbra (CMUC), funded by the European Regional Development Fund through the pro-
gram COMPETE and by the Portuguese Government through the FCT - Fundac¸a˜o para a
Cieˆncia e a Tecnologia under the project PEst-C/MAT/UI0324/2011, and by FCT under
the project PTDC/SAU-ENB/119132/2010) and COMPETE (QREN -FCOMP-01-0124-
FEDER-021014) program.
c CMMSE ISBN: 978-84-616-2723-3Page 127 of 1797
Numerical solution of time-dependent Maxwell’s equations
References
[1] R. Bernardes, T. Santos and J. Cunha-Vaz, Evaluation of Blood-Retinal Bar-
rier Function from Fourier Domain High-Definition Optical Coherence Tomography, in
IFMBE Proceedings 25/XI, 2009, 316–319.
[2] R. Bernardes, T. Santos, P. Serranho, C. Lobo, and J. Cunha-Vaz, Nonin-
vasive evaluation of retinal leakage using OCT, Ophtalmologica 226(2) (2011) 29–36.
[3] B. Bouma and G. Tearney, Handbook of optical coherence tomography, Marcel
Dekker, New York, 2002.
[4] S. Descombes, S. Lanteri and L. Moya, Locally implicit time integration strategies
in a discontinuous Galerkin method for Maxwell’s equations, J. Sci. Comput. 56 (2013)
190–218.
[5] J. Hesthaven and T. Warburton, Nodal High-Order Methods on Unstructured
Grids: I. Time-Domain Solution of Maxwell’s Equations, J. Comput. Phys. 181 (2002)
186–221.
[6] J. S. Hesthaven and T. Warburton, Nodal Discontinuous Galerkin Methods: Al-
gorithms, Analysis, and Applications, Springer Verlag, New York, 2008.
[7] L. Junqueira and J. Carneiro, Basic Histology: Text & Atlas (Junqueira’s Basic
Histology), McGraw-Hill Medical, 2005.
[8] M. Ko¨nig, K. Busch and J. Niegemann, The Discontinuous Galerkin Time-
Domain method for Maxwell’s equations with anisotropic materials, Phot. Nano. Fund.
Appl. 8 (2010) 303–309.
[9] C. Mu¨ller, Foundations of the Mathematical Theory of Electromagnetic Waves,
Springer Verlag, Berlim, 1969.
[10] J. Niegemann, R. Diehl and K. Busch, E cient low-storage Runge-Kutta schemes
with optimized stability regions, J. Comput. Phys. 231 (2012) 364–372.
[11] W. Reed, T. Hill, Triangular mesh methods for the neutron transport equation, Los
Alamos Scientific Laboratory. LA-UR-73-479 (1973).
[12] R.N. Rieben, G.H. Rodrigue and D.A. White, A high order mixed vector finite
element method for solving the time dependent Maxwell equations on unstructured grids,
J. Comput. Phys. 204 (2005) 490–519.
c CMMSE ISBN: 978-84-616-2723-3Page 128 of 1797
A. Arau´jo, S. Barbeiro, L. Pinto, R. Bernardes et al.
[13] Z. S. Sacks, D.M. Kingsland, R. Lee and J. -F. Lee, A perfectly matched
anisotropic absorber for use as an absorbing boundary condition, IEEE Trans. An-
tennas. Propag. 43 (1995) 1460–1463.
[14] J. M. Schmitt, A. Knu¨ttel, and R. F. Bonner, Measurement of optical properties
of biological tissues by low-coherence reflectometry, Appl. Opt. 32 (1993) 6032–6042.
[15] P. Serranho, M. Morgado and R. Bernardes Optical Coherence Tomography: a
concept review In: Optical Coherence Tomography: A Clinical and Technical Update.
R. Bernardes & J. Cunha-Vaz Eds., Springer-Verlag, 2012, 139–156
[16] A. Taflove and S.C. Hagness, Computational Electrodynamic: The Finite-
Di↵erence Time-Domain Method (2nd ed.), Artech House, Norwood, MA, 2000.
[17] K. S. Yee, Numerical solution of initial boundary value problems involving maxwell’s
equations in isotropic media, IEEE Trans. Antennas. Propag. 14 (1966) 302–307.
c CMMSE ISBN: 978-84-616-2723-3Page 129 of 1797


Numerical solution of time-dependent Maxwell’s equations for modeling scattered electromagnetic wave’s propagation

http://estudogeral.sib.uc.pt/jspui/bitstream/10316/44986/1/2013ABPCCMSSB_p.pdf

Numerical solution of time-dependent Maxwell’s equations for modeling scattered electromagnetic wave’s propagation

Abstract

Similar works

Full text

Available Versions

Estudo Geral