Fast Finite Volume Simulation of 3D Electromagnetic Problems with Highly Discontinuous Coefficients.
We consider solving three-dimensional electromagnetic problems in parameter regimes where the quasi-static approximation applies and the permeability, permittivity, and conductivity may vary significantly. The difficulties encountered include handling solution discontinuities across interfaces and accelerating convergence of traditional iterative methods for the solution of the linear systems of algebraic equations that arise when discretizing Maxwell's equations in the frequency domain.The present article extends methods we proposed earlier for constant permeability [E. Haber, U. Ascher, D. Aruliah, and D. Oldenburg, J. Comput. Phys., 163 (2000), pp. 150--171; D. Aruliah, U. Ascher, E. Haber, and D. Oldenburg, Math. Models Methods Appl. Sci., to appear.] to handle also problems in which the permeability is variable and may contain significant jump discontinuities. In order to address the problem of slow convergence we reformulate Maxwell's equations in terms of potentials, applying a Helmholtz decomposition to either the electric field or the magnetic field. The null space of the curl operators can then be annihilated by adding a stabilizing term, using a gauge condition, and thus obtaining a strongly elliptic differential operator. A staggered grid finite volume discretization is subsequently applied to the reformulated PDE system. This scheme works well for sources of various types, even in the presence of strong material discontinuities in both conductivity and permeability. The resulting discrete system is amenable to fast convergence of ILU-preconditioned Krylov methods.We test our method using several numerical examples and demonstrate its robust efficiency. We also compare it to the classical Yee method using similar iterative techniques for the resulting algebraic system, and we show that our method is significantly faster, especially for electric sources.
Introduction
The need for calculating fast, accurate solutions of three-dimensional electromagnetic
equations arises in many important application areas including,
among others, geophysical surveys and medical imaging [29, 32, 2]. Conse-
quently, a lot of effort has recently been invested in finding appropriate numerical
algorithms. However, while it is widely agreed that electromagnetic
phenomena are generally governed by Maxwell's equations, the choice of numerical
techniques to solve these equations depends on parameter ranges and
various other restrictive assumptions, and as such is to a significant degree
application-dependent [20, 32, 2].
The present article is motivated by remote sensing inverse problems, e.g.
in geophysics, where one seeks to recover material properties - especially conductivity
- in an isotropic but heterogeneous body, based on measurements
of electric and magnetic fields on or near the earth's surface. The forward
model, on which we concentrate here, consists of Maxwell's equations in the
frequency domain over a frequency range which excludes high frequencies.
Assuming a time-dependence e \Gamma-!t , these equations are written as
r
r \Theta H \Gamma b
r
r
where - is the magnetic permeability, oe is the conductivity, ffl is the electrical
permittivity,
J s is a known source current density, and ae is the (unknown) volume density
of free charges. In our work we assume that the physical properties - ? 0,
can vary with position, and -ffl! 2 L 2 - 1, where L is a
typical length scale. The electric field E and the magnetic field H are the
unknowns, with the charge density defined by (1c). Note that as long as
redundant and can be viewed as an invariant of the system,
obtained by taking the r\Delta of (1a). The system (1) is defined over a three-dimensional
spatial
domain\Omega\Gamma In principle, the
domain\Omega is unbounded (i.e.
practice, a bounded subdomain of IR 3 is used for numerical
approximations. In this paper we have used the boundary conditions
H \Theta n
although other boundary conditions are possible.
A number of difficulties arise when attempting to find numerical solutions
for this three-dimensional PDE system. These difficulties include handling
regions of (almost) vanishing conductivity, handling different resolutions in
different parts of the spatial domain, handling the multiple scale lengths over
which the physical properties can vary, and handling regions of highly varying
conductivity, magnetic permeability, or electrical permittivity where jumps
in solution properties across interfaces may occur.
On the other hand, the nature of the data (e.g. measurements of the
electric and/or magnetic fields at the surface of the earth) is such that one
cannot hope to recover to a very fine detail the structure of the conductivity oe
or the permeability -. We therefore envision, in accordance with the inverse
problem of interest, a possibly nonuniform tensor product grid covering the
domain\Omega\Gamma where b oe and - are assumed to be smooth or even constant inside
each grid cell, but they may have significant jump discontinuities that can
occur anywhere
in\Omega across cell interfaces. The source J s is also assumed not
to have jumps across interfaces. The relative geometric simplicity resulting
from this modeling assumption is key in obtaining highly efficient solvers for
the forward problem.
Denoting quantities on different sides of an interface by subscripts 1 and
2, it can be easily shown [33] that across an interface
n \Theta
These conditions imply that neither E nor H are continuous in the normal
direction when b oe, resp. -, have a jump discontinuity across a cell face,
and likewise, b oeE and -H are not necessarily continuous in tangential di-
rections. Care must therefore be exercised when numerical methods are employed
which utilize these variables if they are to be defined where they are
double-valued.
By far the most popular discretization for Maxwell's equations is Yee's
method [37] (see discussions and extensions in [32, 26, 19]). This method
employs a staggered grid, necessitating only short, centered first differences
to discretize (1a) and (1b). In the more usual application of this method, the
electric field components are envisioned on the cell's edges and the magnetic
field components are on the cell's faces - see Fig. 1. It is further possible to
H z
x
y
z

Figure

1: A staggered discretization of E and H in three dimensions: E-
components on edges, H-components on faces.
eliminate the components of the magnetic field from the discrete equations,
obtaining a staggered discretization for the second order PDE in E,
r \Theta (- \Gamma1 r \Theta E) \Gamma
Related methods include the finite integration technique and certain mixed
finite element methods [35, 5, 25]. Although these methods are often presented
in the context of time-domain Maxwell's equations the issues arising
when applying an implicit time-discretization (a suitable technique under our
model assumptions) are often similar to the ones we are faced with here.
The popularity of Yee's method is due in part to its conservation properties
and other ways in which the discrete system mimics the continuous
system [19, 15, 5, 4]. However, iterative methods to solve the discrete system
may converge slowly in low frequencies, due to the presence of the rich,
nontrivial null space of the curl operator, and additional difficulties arise
when highly discontinuous coefficients are present [29, 24, 16]. There are two
major reasons for these difficulties. First, the conductivity can essentially
vanish (for example, in the air, which forms part of \Omega\Gamma/ from an analytic
perspective, the specific subset of Maxwell's equations used typically forms
an almost-singular system in regions of almost-vanishing b oe. Even in regions
where the conductivity is not close to vanishing, the resulting differential
operator is strongly coupled and not strongly elliptic [6, 1]. Second, in cases
of large jump discontinuities, care must be taken to handle H and b
carefully, since these are located as in Fig. 1 where they are potentially discontinuous

In [1], we addressed the often slow convergence of iterative methods when
used for the equations resulting from the discretization of (5) by applying
a Helmholtz decomposition first, obtaining a potential formulation with a
Coulomb gauge condition. This change of variables (used also in [3, 12,
22, 28], among many others) splits the electric field into components in the
active and the null spaces of the curl operator. A further reformulation,
reminiscent of the pressure-Poisson equation for the incompressible Navier-Stokes
equations [14, 31], yields a system of strongly elliptic, weakly coupled
PDEs, for which more standard preconditioned Krylov space methods are
directly applicable.
In [15], we further addressed possible significant jumps in the conductivity
while - is assumed constant, by employing a finite volume discretization on
a staggered grid, akin to Yee's method with the locations of E- and H-
components exchanged, as in Fig. 2. The normal components of E are now
double-valued, but this is taken care of in an elegant way by the Helmholtz
decomposition of E and by introducing the (generalized) current
x
y
z

Figure

2: A staggered discretization of E and H in three dimensions: E-
components on faces, H-components on edges.
into the equations. The curl operators in (5) are replaced by the vector
Laplacian according to the vector identity
r \Theta r\Theta
for sufficiently smooth vector functions (not E).
In this paper we generalize our approach from [1, 15] to the case where
the magnetic permeability - may be highly discontinuous as well. This is
a realistic case of interest in geophysical applications, although usually the
jump in conductivity dominates the jump in permeability. Now the roles
of E and H are essentially dual, and it is possible to apply a Helmholtz
decomposition to either E or H, keeping the other unknown vector function
intact. We choose to decompose the electric field E, referring to Fig. 2 for
the locations of the H-unknowns in the ensuing discretization. The major
departure from our previous work is in the fact that the identity (7) does not
directly extend for the operator r \Theta (- \Gamma1 r\Theta ) appearing in (5). We can,
however, stabilize this operator by subtracting r(-
this forms the basis for our proposed method. In cases of constant magnetic
permeability or electric conductivity the formulation can be reduced to our
previous formulation in [15] or a variant thereof.
Our approach to dealing with possible discontinuities can be viewed as
using fluxes (which are continuous on cell faces) and vorticities (which are
continuous on cell edges). The introduction of such unknowns is strongly
connected to mixed finite elements which are used for highly discontinuous
problems [7, 8, 16].
The paper is laid out as follows. In Section 2, we reformulate Maxwell's
equations in a way which enables us to extend our methods. The resulting
system is amenable to discretization using a finite-volume technique described
in Section 3.
The extension and generalization of our approach from [1] through [15] to
the present article is not without a price. This price is an added complication
in the sparsity structure of the resulting discrete system and a corresponding
added cost in the iterative solution of such systems. We briefly describe
the application of Krylov space methods to solve the system of algebraic
equations in Section 4. We use BICGSTAB (see, e.g., [30]) together with one
of two preconditioners: an incomplete block LU-decomposition (which is a
powerful preconditioner in the case of diagonally dominant linear systems)
and SSOR. The system's diagonal dominance is a direct consequence of our
analytic formulation.
We present the results of numerical experiments in Section 5 and compare
results obtained using our method with those obtained using a more traditional
Yee discretization. If the source is not divergence-free, as is the case
for electric (but not magnetic) sources, then our method is better by more
than two orders of magnitude. The method works well also for a case where
the problem coefficients vary rapidly. We conclude with a short summary
and further remarks.
2 Reformulation of Maxwell's Equations
Maxwell's equations (1a) and (1b) can be viewed as flux balance equations,
i.e. each term describes the flux which arises from a different physical con-
sideration, and the equations are driven by the conservation of fluxes. (In
fact, this was how they were originally developed [23].) Therefore, in both
(1a) and (1b) we have a flux term which should be balanced. In (1b) the
generalized current density b
J defined in (6) is balanced with the source and
the flux which arise from magnetic fields, and in (1a) the magnetic flux
is balanced with the flux which arises from electric fields. In our context
these fluxes are well-defined on cell faces but they may be multi-valued at
cell edges. 1
Furthermore, the leading differential operator in (5), say, has a nontrivial
null space. Rather than devising iterative methods which directly take
this into account (as, e.g., in [2, 16]), we transform the equations before
discretization.
We decompose E into its components in the active space and in the null
space of the curl operator:
r
We could decompose H instead in a similar way, but we would not decompose
both. Here we have chosen to concentrate on the decomposition of E.
Substituting equations we obtain
r \Theta A \Gamma
r \Theta H \Gamma b
r
Furthermore, (5) becomes
r \Theta (- \Gamma1 r \Theta
r
Note that across an interface between distinct conducting media we have,
in addition to (4),
n \Theta
ae s in (12c) is an electric surface charge density. These conditions and
the differential equations (1) imply that while b
H(curl;\Omega\Gamma4 see,
e.g., [13].
not. Moreover, rOE inherits the discontinuity of E \Delta n, while A is continuous,
and both r \Delta A and r \Theta A are bounded (cf. [13]).
In [15] we had the relation r \Theta r \Theta A = \Gammar 2 A holding. However,
when - varies, the identity (7) does not extend directly, and we must deal
with the null space of r \Theta A in a different way.
Let us define the Sobolev spaces
equipped with the usual norm (see, e.g., [13])
kvk
(the
(\Omega\Gamma4 and [L
are used on the right hand side of (13b)),
and
Green's formula yields, for any u 2
(r \Theta (-
where the usual notation for inner product in L
2(\Omega\Gamma and [L
2(\Omega\Gamma2 3 is used.
Thus, for any u 2 W
(r \Theta (- \Gamma1 r \Theta u); u)
We may therefore stabilize (11a) by subtracting a vanishing term:
r \Theta (- \Gamma1 r \Theta
obtaining a strongly elliptic operator for A, provided A 2 W
0(\Omega\Gamma4 The
latter condition is guaranteed by the choice (17b) below.
A similar stabilization (or penalty) approach was studied with mixed
success in [10, 20]. However, our experience and thus our recommendation
are more positive because of the discretization we utilize. We elaborate
further upon this point in the next section.
Using (10c), we can write (10b) as
r \Theta H \Gamma
This may be advantageous in the case of discontinuous -, similarly to the
mixed formulation used for the simple div-grad system
r
in [7, 27, 8, 11] and elsewhere.
Our final step in the reformulation is to replace, as in [15], the gauge
condition (10c) on A by an indirect one, obtained upon taking r\Delta of (15a)
and simplifying using (15b) and (10c). This achieves only a weak coupling
in the resulting PDE system. We note that replacing of the gauge condition
(10c) is similar to the pressure Poison equation in CFD [31, 14]. The complete
system, to be discretized in the next section, can now be written as
r \Theta H \Gamma
In order to complete the specification of this system we must add appropriate
Boundary Conditions (BCs). First, we note that the original BC (3)
can be written as
(r \Theta A) \Theta n
An additional BC on the normal components of A is required for the
Helmholtz decomposition (9) to be unique. Here we choose (corresponding
to (13c))
A
This, together with (9), determines A for a given E.
Moreover, since (16d) was obtained by taking the r\Delta of (16a) additional
are required on either @OE=@n or OE. For this we note that the original
together with the PDE (1b) imply also
The latter relation (17c), together with (17b), implies @OE
@n
0 at the bound-
ary. 23
The above conditions determine OE up to a constant. We pin this constant
arbitrarily, e.g. by requiring
Z
Finally, we note that (17c) together with (16a) and (3), imply in turn
that @/
@n
0 at the boundary. Since (16a) and (16d) imply that
we obtain / j 0, and thus retrieve (10c), by pinning / down to 0 at one
additional point.
The system (16) subject to the boundary conditions (17) and / pinned
at one point is now well-posed.
3 Deriving a discretization
As in [15], we employ a finite volume technique on a staggered grid, where b
and A are located at the cell's faces, H is at the cell's edges, and OE and / are
located at the cell's center. Correspondingly, the discretization of (16d) and
(16c) are centered at cell centers, those of (16e) and (16a) are centered at cell
faces, and that of (16b) is centered at cell edges. The variables distribution
over the grid cell is summarized in Table 1.
To approximate r \Delta u for
integrate at first over the cell using the Gauss divergence theoremje i;j;k j
Z
e i;j;k
r \Delta u
Z
@e i;j;k
In cases where the original BC is different from (3) we still use the boundary condition
as an asymptotic value for an infinite domain. Alternatively, note the
possibility of applying the Helmholtz decomposition to H , although generally there are
good practical reasons to prefer the decomposition (9) of E.
3 In our geophysical applications J s \Delta n vanishes at the boundary.
A x
A y
A z
J z
J z
H z

Table

1: Summary of the discrete grid functions. Each scalar field is approximated
by the grid functions at points slightly staggered in each cell e i;j;k of
the grid.
and then use midpoint quadrature on each face to evaluate each component
of the surface integrals appearing on the right-hand side above. Thus, define
and express the discretization of (16d) and (16c) on each grid cell as
Note that we are not assuming a uniform grid:
each cell may have different widths in each direction. The boundary conditions
are used at the end points of (19).
Next, consider the discretization at cell faces. Following [15], we define
the harmonic average of boe between neighboring cells in the x-direction by
where h x
)=2. If b oe is assumed constant over each
cell, this integral evaluates to
2boe i;j;k
2boe i+1;j;k
Then, the resulting approximation for the x-component of (16e) is [15]
A x
Next, we discretize (16a) as in [37]. Writing the x-component of these
equations,
where we denote J discretization centered at the center of
the cell's x-face results:
H z
Similar expressions to those in (20c) and (21) can be derived in the yand
z-directions. The boundary conditions (3) are used to close (21). Using
(20c) we can eliminate b
J from (19a) and obtain a discrete equation in which
the dominant terms all involve OE. The resulting stencil for OE has 7 points.
We also apply the obvious quadrature for the single condition (17d).
Finally we discretize the edge-centered (16b). Consider, say, the x-component
of (16b), written as
Integrating this equation over the surface of a rectangle with corners
the expression on the left hand side is integrated using the Gauss curl theo-
rem, and - on the right hand side is averaged in both directions to obtain a
value on the edge. We do not divide through by - before this integration because
we wish to integrate the magnetic flux, which is potentially less smooth
around an edge than the magnetic field. This yields
Z z k+1
z k
If - is assumed to be constant over each cell, this integral evaluates to
[h y
h z
h z
h z
h z
Then, the resulting approximation for the x-component of (16b) is
iA y
h z
A z
(22c)
Using (19b) as well as (22c) and similar expressions derived in the y- and
z-directions, we substitute for H and / in (21) and obtain a discrete system
of equations for A. The resulting stencil for A has 19 points and the same
structure as for the discretization of the operator
r \Theta (- \Gamma1 r\Theta
The difference between this discretization and a direct discretization of the
latter is that - at the interface is naturally defined as an arithmetic average
and not a harmonic average.
The discretization described above can be viewed as a careful extension
of Yee's method, suitable for discontinuous coefficients and amenable to fast
iterative solution methods. It is centered, conservative and second order accu-
rate. Note that throughout we have used a consistent, compact discretization
of the operators r\Delta ; r\Theta and r. We can denote the corresponding discrete
operators by r h and r h and immediately obtain the following identities
(cf. [18, 17]),
These are of course analogues of vector calculus identities which hold for
sufficiently differentiable vector functions. The boundary conditions (17) are
discretized using these discrete operators as well.
Next, note that upon applying r h \Delta to (21) and using (23b) and (19a) we
obtain
Moreover, from (21) and (17c), the discrete normal derivative of / vanishes
at the boundaries as well. Setting / to 0 arbitrarily at one point,
then determines that / j 0 throughout the domain (as a grid function).
We obtain another conservation property of our scheme, namely, a discrete
divergence-free A:
Recall the stabilizing term added in (14). For the exact solution this term
obviously vanishes. Now, (24) assures us that the corresponding discretized
term vanishes as well (justifying use of the term 'stabilization', rather than
'penalty'). This is not the case for the nodal finite element method which was
considered in [20, 10]. For an approximate solution that does not satisfy (24)
the stabilization term may grow in size when - varies over a few orders of
magnitude, or else r h \Delta A grows undesirably in an attempt to keep -
approximately constant across an interface [20].
Our particular way of averaging across discontinuities, namely, arithmetic
averaging of - at cell edges and harmonic averaging of b
oe at cell faces, can
be important. The averaging can be viewed as a careful approximation of
the constitutive relationship for discontinuous coefficients. To show that, we
look first at the relation
across a face whose normal direction is x. This flux flows in series, and
therefore an approximate b oe that represents the bulk property of the flow
through the volume is given by the harmonic average (corresponding to an
arithmetic average of the resistivities). Next, we look at the relation
where - is an edge variable and B x is the flux through four cells which
share that edge. Here the flow is in parallel, which implies that we need to
approximate - on the edge by an arithmetic average.
Note also that if we use the more common implementation of Yee's
method (i.e. H on the cell's faces and E on the cell's edges) then the
roles of b oe and - interchange and we need to average - harmonically and b oe
arithmetically.
4 Solution of the discrete system
After the elimination of H , / and b
J from the discrete equations we obtain
a large, sparse system for the grid functions corresponding to A and OE:
A
OE
b A
Here
is the result of the discretization of the operator r \Theta -
corresponds to the discretization of the operator r\Theta ; D likewise corresponds
to the discretization of the operator the diagonal matrix S results from
the discretization of the operator b oe(\Delta); M and M c similarly arise from the
at cell edges and at cell centers respectively; and H
represents the discretization of r \Delta (boer(\Delta)). In regions of constant - the
simplifies into a discretization of the vector Laplacian. The blocks
are weakly coupled through interfaces in -. A typical sparsity structure of
H 1 for variable -, as contrasted with constant -, is displayed in Fig. 3. The
structure of the obtained system is similar to that in [15], although the main
block diagonal is somewhat less pleasant.
Note that as long as the frequency is low enough that the diffusion number
satisfies
(where h is the maximum grid spacing) the matrix is dominated by the
diagonal blocks. This allows us to develop a block preconditioner good for
low frequencies ! based on the truncated ILU (ILU(t)) decomposition of the
major blocks [30]. Thus, we approximate K by the block diagonal matrix
and then use ILU(t) to obtain a sparse approximate decomposition of the
K. This decomposition is used as a preconditioner for the Krylov
solver. Note that although K is complex the approximation -
K is real and
therefore we need to apply the ILU decomposition only to two real matrices,
which saves memory.
Note that the block approximation (26) makes sense for our discretization
but not for the direct staggered-grid discretization of (1). Thus, the
Sparsity structure - constant -
Sparsity structure - variable -

Figure

3: Sparsity structure of the matrix H 1 corresponding to variable -
(right) and to constant - (left).
reformulation and subsequent discretization of Maxwell's equations allow us
to easily obtain a good block preconditioner in a modular manner, while the
discretization of (1) does not.
5 Numerical examples
Our goal in this section is to demonstrate the efficacy of our proposed method
and to compare it to the more common discretization of the system (1), utilized
e.g. in [24, 29], using standard Krylov-type methods and preconditioners
for solving the resulting discrete systems. We vary the type of source
used, the size of jumps in the coefficients oe and -, the preconditioner and
the grid size.
In the tables below, 'iterations' denotes the number of BICGSTAB iterations
needed for achieving a relative accuracy of
the number of giga-flops required; the SSOR parameter (when used) equals
and the ILU threshold (when used) equals 10 \Gamma2 . The latter threshold is
such that in our current Matlab implementation iterations involving these
two preconditioners cost roughly the same.
5.1 Example Set 1
We derive the following set of experiments. Let the air's permeability be
its permittivity ffl We assume a
cube of constant conductivity oe c and permeability - c embedded in an otherwise
homogeneous earth with conductivity oe permeability
Fig. 4. In a typical geophysical application, the conductivity
may range over four orders of magnitude and more, whereas the permeability
rarely changes over more than one order of magnitude. Therefore, we
experiment with conductivity oe c ranging from 10 \Gamma2 S/m to 10 3 S/m and
permeability - c ranging from - 0 to 1000- 0 .
-0.4

Figure

4: The setting of our first numerical experiment. A cube of conductivity
oe c and permeability - c is embedded inside a homogeneous earth with
conductivity oe e and permeability - e . Also, in the air
We experiment with two different sources: (i) a magnetic source (a plane
wave); and (ii) an electric dipole source in the x direction centered at (0; 0; 0).
The fact that the first source is magnetic implies that it is divergence-free.
This source lies entirely in the active space of the curl operator. In contrast,
the electric source is not divergence-free.
Both sources are assumed to oscillate with different frequencies ranging
from 1 to 10 6 Hz. The solution is obtained, unless otherwise noted, on a
nonuniform tensor grid (see Fig. 4) consisting of cells. There are 95232
unknowns corresponding to this grid and 128000 (complex)
unknowns. We then solve the system using the method described in the
previous section.
5.1.1 Example 1a.
In order to be able to compare the resulting linear algebra solver with that
corresponding to Yee's method we discretize the system
r \Theta [- \Gamma1 r \Theta ( b
using the staggered grid depicted in Fig. 2, i.e., where b
J is on cells' faces,
which is similar to the discretization in [24]. This yields the discrete system
for the unknowns vector e corresponding to grid values of b
J=boe, where the
matrices C; M and S are defined in Section 4 and - b depends on the source.
In order to solve this system as well as ours we use BICGSTAB and an
preconditioner. The comparison between the methods for the case
different frequencies is summarized in Table 2.
Electric source Magnetic source
iterations operations iterations operations
6.3 728 99 113 6.4 7.3

Table

2: Iteration counts and computational effort for our method (A; OE)
and the traditional implementation of Yee's method (applied to E, or b
J
Example Set 1 using both an electric source and a magnetic source.

Table

2 shows that our method converges in a moderate number of iterations
for both sources, despite the presence of significant jumps in - and
oe. On the other hand, the more traditional discretization performs poorly
for the electric source and reasonably well for the magnetic source. Slow
convergence of the direct staggered discretization of Maxwell's equations in
the case of an electric source was also reported in [29], where E was defined
on the grid's edges.
These results clearly show an advantage of our formulation over the original
Yee formulation, even for a simple preconditioning, especially for electric
sources and in low frequencies. In such circumstances, the discretized first
term on the left hand side of (27) strongly dominates the other term, and
the residual in the course of the BICGSTAB iteration has a nontrivial component
in the null space of that operator; hence its reduction is very slow.
The magnetic source, on the other hand, yields a special situation where, so
long as the discrete divergence and the roundoff error are relatively small,
the residual component in the null space of the leading term operator is also
small, hence the number of iterations using the traditional method is not
much larger than when using our method.
5.1.2 Example 1b.
Next, we test the effect of discontinuities on our method. We use the electric
source and record the number of BICGSTAB iterations needed for our
method to converge for various values of ~
using block-ILU(t) preconditioning as described in the previous section. The
results are summarized in Table 3.
Note that large jump discontinuities in oe do not significantly affect the
rate of convergence of the iterative linear system solver for our method, but
large jump discontinuities in - have a decisive effect. Results in a similar
spirit were reported in [16] regarding the effect of discontinuities in - on a
specialized multigrid method for an edge-element discretization. However,
even for large discontinuities in - the number of iterations reported in Table
3 remains relatively small compared with similar experiments reported in
[29, 9]. We attribute the increase in the number of iterations as the jump in
- increases in size to a corresponding degradation in the condition number
of K in (25). This degradation, however, does not depend strongly on grid
size, as we verify next.
Table

3: Iteration counts for different frequencies, conductivities and permeabilities
in Example Set 1. The conductivity/permeability structure is a
cube in a half-space.
5.1.3 Example 1c.
In the next experiment, we use the cube model with the electric source to
evaluate the influence of the grid on the number of iterations. We fix
and test our method on a case with a modest coefficient jump ~
and on a
case with a large jump ~
set of uniform grids in the interval [\Gamma1; 1] 3
is considered. For each grid we record the resulting number of iterations
using both the SSOR and the block ILU preconditioners. The results of this
experiment are gathered in Table 4.
We observe that the number of iterations increases as the number of un-
Grid size ILU SSOR ILU SSOR

Table

4: Iteration counts for different grids, for two sets of problem coefficients
and using two preconditioners.
knowns increases. The increase appears to be roughly proportional to the
number of unknowns to the power 1=3. The growth in number of iterations
as a function of grid size is also roughly similar for both preconditioners,
although the block ILU requires fewer iterations (about 1=4 for ~
for each grid. However, ILU requires more memory than SSOR, which may
prohibit its use for very large problems. The increase rate is also similar, as
expected, for both values of ~
-. Thus, the increase in number of iterations as
a function of ~
essentially does not depend on the grid size. Practically, how-
ever, this increase is substantial and may be hard to cope with for (perhaps
unrealistically) large values of ~
- using the present techniques.
5.2 Example 2
In our next experiment we consider a more complicated earth structure. We
employ a random model, which is used by practitioners in order to simulate
stochastic earth models [21]. Two distinct value sets oe
probability P are assumed for the conductivities and permeabilities: for each
cell, the probability of having values oe and the probability of having
values . This can be a particularly difficult model to work with,
as the conductivity and permeability may jump anywhere from one cell to
the next, not necessarily just across a well-defined manifold. A cross-section
of such a model is plotted in Fig. 5. We then carry out experiments as before
for frequencies ranging from 0 to 10 6 Hz.
We use the random model with different conductivities, permeabilities
and frequencies and the electric source, with on the domain [\Gamma1; 1] 3
(in km). We employ a uniform grid of size 44 3 and use both the block
ILU and the SSOR preconditioners. The results of this experiment in terms
Conductivity slice log(S/m)

Figure

5: The setting of Example Set 2.
of iteration counts are summarized in Table 5. The results show that our
solution method is effective even for highly varying conductivity. As before,
the method deteriorates when very large variations in - are present.
We can also see from Table 5 that the block ILU preconditioner works
very well for low frequencies, but it is not very effective for high frequencies.
It is easy to check that in all cases where the block ILU preconditioner fails
to achieve convergence (denoted 'nc', meaning failure to achieve a residual of
\Gamma7 in 800 iterations) the maximum grid spacing h satisfies !-oe h AE 1. In
such a case the discretization of the leading terms of the differential operator
no longer yields the dominant blocks in the matrix equations (25), and therefore
our block ILU preconditioner fails. Thus, for high frequency and high
conductivity we require more grid points in order for this preconditioner to
be effective. This is also consistent with the physics, as the skin depth [36]
decreases and the attenuated wave can be simulated with fidelity only on a
finer grid.
6 Summary and further remarks
In this paper we have developed a fast finite volume algorithm for the solution
of Maxwell's equations with discontinuous conductivity and permeability.
The major components of our approach are as follows.
ffl Reformulation of Maxwell's equations: The Helmholtz decomposition is
applied to E; then a stabilizing term is added, resulting in a strongly elliptic
system; the system is written in first order form to allow flexibility
ILU SSOR ILU SSOR ILU SSOR

Table

5: Iteration counts for different frequencies, conductivities and permeabilities
for Example Set 2.
in the ensuing discretization; and finally, the divergence-free Coulomb
gauge condition is eliminated using differentiation and substitution,
which yields a weakly coupled PDE system enabling an efficient preconditioner
for the large, sparse algebraic system which results from
the discretization.
ffl Discretization using staggered grids, respecting continuity conditions
and carefully averaging material properties across discontinuities. For
this discretization, the stabilizing term vanishes at the exact solution of
the discrete equations, which is important for cases with large contrasts
in -.
ffl Solution of the resulting linear system using fast preconditioned Krylov
methods.
The resulting algorithm was tested on a variety of problems. We have shown
dramatic improvement over the more standard algorithm when the source is
electric. Good performance was obtained even when the coefficients - and
oe were allowed to vary rapidly on the grid scale - a case which should prove
challenging for multigrid methods.
The project that has motivated us is electromagnetic inverse problems
in geophysical prospecting [34]. Solving the forward problem, i.e. Maxwell's
equations as in the present article, is a major bottleneck for the data inversion
methods - indeed it may have to be carried out dozens, if not hundreds, of
times for each data inversion. Thus, extremely fast solvers of the problem
discussed in our paper are needed. Based on the algorithm described here an
implementation has been carried out which solves realistic instances of this
forward problem in less than two minutes on a single-processor PC, enabling
derivation of realistic algorithms at low cost for the inverse problem.

Acknowledgment

We wish to thank Drs. Doug Oldenburg and Dave
Moulton for fruitful discussions and an anonymous referee for valuable comments
on the first two sections of our exposition.



--R

A method for the forward modelling of 3D electromagnetic quasi-static problems
Adaptive multilevel methods for edge element discretizations of maxwell's equations.

Whitney forms: a class of finite elements for three-dimensional computations in electromagnetism

Multigrid solution to elliptic flow problems.
Mixed and hybrid finite element methods.
Finite Elements
Numerical algorithms for the FDITD and FDFD simulation of slowly varying electromagnetic feilds.
Theoretical and numerical difficulties in 3D vector potential methods in finite element magnetostatic compu- tations
Black box multigrid.
Geomagnetic induction in a heterogenous sphere: Azimuthally symmetric test computations and the response of an undulating 660-km discontinuity
Finite Element Methods for Navier-Stokes Equations
On pressure boundary conditions for the incompressible Navier-Stokes equations
Fast modelling of 3D electromagnetic using potentials.
Multigrid method for maxwell's equations.
Natural discretizations for the divergence gradient and curl on logically rectangular grids.
The orthogonal decomposition theorems for mimetic finite difference methods.
Mimetic discretizations for Maxwell's equations and equations of magnetic diffusion.
The Finite Element Method in Electromagnetics.
Electrical conduction in sandstones: Examples of two dimensional and three dimensional percolation.
A scalar-vector potential solution for 3D EM finite-difference modeling
Theoretical concepts in physics.
Three dimensional magnetotelluric modeling using difference equations - theory and comparison to integral equation solutions
An analysis of finite volume
A convergence analysis of Yee's scheme on nonuniform grids.
Numerical implementation of nodal methods.
Solutions of 3D eddy current problems by finite elements.

Iterartive Methods for Sparse Linear Systems.
A multigrid solver for the steady state navier-stokes equations using the pressure-poisson formulation
Computational Electrodynamics: the Finite-Difference Time-Domain Method
Diadic Green Functions in Electromagnetics.
Electromagnetic theory for geophysical applications.
Time domain electromagnetic field computation with finite difference methods.
Inversion of Magnetotelluric Data for a One Dimensional Conductivity
Numerical solution of initial boundary value problems involving maxwell's equations in isotropic media.
--TR

--CTR
Eldad Haber , Stefan Heldmann, An octree multigrid method for quasi-static Maxwell's equations with highly discontinuous coefficients, Journal of Computational Physics, v.223 n.2, p.783-796, May, 2007
krylov methods;mixed methods;coulomb gauge;solution discontinuities;finite volume;helmholtz decomposition;maxwell's equations;preconditioning