University of Maryland Department of Computer Science TR-5040
University of Maryland Institute for Advanced Computer Studies TR-2014-06
July 2014
A STOCHASTIC APPROACH TO UNCERTAINTY IN THE
EQUATIONS OF MHD KINEMATICS∗
EDWARD G. PHILLIPS† AND HOWARD C. ELMAN‡
Abstract. The magnetohydodynamic (MHD) kinematics model describes the electromagnetic
behavior of an electrically conducting fluid when its hydrodynamic properties are assumed to be
known. In particular, the MHD kinematics equations can be used to simulate the magnetic field
induced by a given velocity field. While prescribing the velocity field leads to a simpler model than
the fully coupled MHD system, this may introduce some epistemic uncertainty into the model. If the
velocity of a physical system is not known with certainty, the magnetic field obtained from the model
may not be reflective of the magnetic field seen in experiments. Additionally, uncertainty in physical
parameters such as the magnetic resistivity may affect the reliability of predictions obtained from
this model. By modeling the velocity and the resistivity as random variables in the MHD kinematics
model, we seek to quantify the effects of uncertainty in these fields on the induced magnetic field.
We develop stochastic expressions for these quantities and investigate their impact within a finite
element discretization of the kinematics equations. We obtain mean and variance data through
Monte-Carlo simulation for several test problems. Toward this end, we develop and test an efficient
block preconditioner for the linear systems arising from the discretized equations.
Key words. magnetohydrodynamics, uncertainty quantification, kinematics equations, iterative
methods
1. Introduction. Magnetohydrodynamics (MHD) is the study of the interaction
between electrically conducting fluids and magnetic fields. The MHD model applies
to a range of fluids including plasmas, liquid metals, and brine. The equations govern-
ing MHD dynamics result from a coupling of the equations of MHD kinematics and
the Navier-Stokes equations for fluid flow through the Lorentz force. In this work,
we consider the kinematics equations, which govern the influence of the fluid flow
on the magnetic field. These equations constitute a component block within fully
coupled MHD simulations [3, 17]. Furthermore, solution of these equations is also
required in operator splitting techniques that alternate between solving the Navier-
Stokes equations and the kinematics equations [16]. The kinematics equations are
also of particular interest in the field of kinematic dynamo theory, in which the ratio
of the Lorentz force to inertia is assumed to be small [12]. In this case, the velocity
can be prescribed, and the generation of the magnetic energy induced by the flow can
be studied. Kinematic simulations can be used to model MHD generators, in which
plasmas act as conductors to generate electric currents, as well as natural dynamos
such as the sun and the geodynamo. They are of primary interest in investigating
whether a given flow profile can sustain dynamo action.
When the velocity field is prescribed, this simplifies the MHD equations, but it
may also introduce some epistemic uncertainty into the model. The flow properties
of the fluid may not be known on the interior of the domain. Additionally, there are
aspects of the physical model that motivate incorporation of small-scale uncertainty.
For instance, the large-scale mean flow of the earth’s outer core cannot account for
the magnitude of the earth’s magnetic field. In geodynamo theory, it is proposed that
∗This work was supported in part by the U.S. Department of Energy under grant DE-SC0009301
and the U.S. National Science Foundation under grant DMS1115317.
†Applied Mathematics & Statistics, and Scientific Computation Program, University of Maryland,
College Park, MD (egphillips@math.umd.edu)
‡Department of Computer Science and Institute for Advanced Computer Studies, University of
Maryland, College Park, MD (elman@cs.umd.edu)
1
small-scale turbulent behavior can give rise to a large-scale magnetic field through the
α-effect [5]. Furthermore, the distribution of material properties may be uncertain
in physical applications. When multiple fluids are present, such as when multiple
liquid metals are mixing together, the magnetic resistivity will not be homogeneous
throughout the domain and may vary over orders of magnitude. Because the resistivity
can have a strong influence on such physical systems, including changing the topology
of the magnetic field, we are interested in how uncertain heterogeneous distributions
of the resistivity may affect the induced magnetic energy.
In this study, we explore these issues by mathematically simulating uncertainty in
both the velocity field and the resistivity within the MHD kinematics model. In this
model, we treat the uncertain quantities as random fields correlated in space. We will
obtain mean and variance data through Monte-Carlo simulation. In addition, because
each Monte-Carlo trial requires the solution of linear systems with randomly varying
dynamics, we develop and explore efficient and robust solvers for discrete kinematics
systems.
Thus, we consider two issues: the impact of uncertainty of velocity and resistivity
on statistical properties of the magnetic fields modeled by the equations of MHD kine-
matics, together with efficient computational algorithms for computing these quan-
tities. The remainder of the paper is structured as follows. In Section 2, we will
derive a finite element formulation for the deterministic equations of MHD kinemat-
ics. Section 3 is devoted to the incorporation of uncertainty into the model. In this
section, we describe a means of modeling the uncertainty in both the resistivity and
the velocity field and apply the model to representative test problems. In Section 4,
we propose, analyze, and test a block preconditioner for solving the linear systems
arising in our model. Finally, we will draw conclusions in Section 5.
2. A Finite Element Formulation. The steady-state kinematics of MHD are
governed by Maxwell’s equations
∇×
(
1
µ
~B
)
= ~j, (2.1a)
∇ · (ε ~E) = ρc, (2.1b)
∇× ~E = ~0, (2.1c)
∇ · ~B = 0, (2.1d)
and Ohm’s law
~j = σ( ~E + ~u× ~B), (2.2)
on a domain D ⊂ Rd, d = 2 or 3 (plus appropriate boundary conditions). The
unknowns here are the magnetic induction ~B, the electric field ~E, and the current
density ~j; the charge density ρc can be regarded as an auxiliary variable obtained after
computing ~E. We will regard the fluid velocity ~u as given. For many applications, the
electric permittivity ε and the magnetic permeability µ do not vary significantly for
different fluids [6], so we let ε and µ be fixed scalar constants over the whole domain.
However, because the heterogeneities of the electric conductivity σ can be large for
different liquid metals, we consider σ to be a prescribed, not necessarily constant
function on D.
2
We consider the boundary conditions
~B × ~n = ~q, (2.3a)
~E · ~n = k, (2.3b)
on ∂D. We choose these conditions over the alternative (prescribing ~B · ~n and ~E ×
~n) because the requirement on the tangential component of ~B is then the natural
Dirichlet condition for the curl-conforming edge elements employed to discretize ~B.
A standard simplification of equations (2.1) and (2.2) is obtained by eliminating
the variables ~j and ~E, yielding the following equations for the kinematics of MHD in
terms of ~B:
∇×
(
η
µ
∇× ~B
)
−∇× (~u× ~B) = ~0, (2.4a)
∇ · ~B = 0, (2.4b)
on D, where η = 1/σ is the magnetic resistivity. A boundary condition such as (2.3a)
is required to complete this system. After ~B is obtained from solving equations (2.4),
~E,~j, and ρc can be recovered.
As stated, the equations (2.4) are over-determined because there are d+ 1 equa-
tions in d unknowns. In order to make the system well-defined without changing the
solution ~B, we introduce a Lagrange multiplier r (we refer to this variable as the
magnetic pseudo-pressure), and consider the equations
∇×
(
η
µ
∇× ~B
)
−∇× (~u× ~B) +∇r = ~0, (2.5a)
∇ · ~B = 0, (2.5b)
with the boundary conditions (2.3a) and r = 0 on ∂D. It can be shown that (2.5)
admits the same solution ~B as (2.4) by taking the divergence of equation (2.5a). This
yields ∆r = 0 on D, which, with the zero Dirichlet condition on r, implies that r = 0
on D.
For developing a weak formulation for this problem, we consider the spaces
V0 = {~C ∈ H(curl,D)|~C × ~n = ~0 on ∂D}, (2.6a)
V~q = {~C ∈ H(curl,D)|~C × ~n = ~q on ∂D}, (2.6b)
Q0 = H
1
0 (D). (2.6c)
Multiplying the equations (2.5) by test functions ~C ∈ V0 and s ∈ Q0, and integrating
by parts, we obtain the following weak formulation: Find ( ~B, r) ∈ V~q ×Q0 such that
a( ~B, ~C) + c( ~B, ~C) + b(~C, r) = 0 (2.7a)
b( ~B, s) = 0 (2.7b)
for all (~C, s) ∈ V0 ×Q0, where the bilinear forms are defined as
a( ~B, ~C) =
(
η
µ
∇× ~B,∇× ~C
)
, (2.8a)
c( ~B, ~C) = −
(
~u× ~B,∇× ~C
)
, (2.8b)
b( ~B, s) = ( ~B,∇s). (2.8c)
3
By integrating by parts the divergence of the magnetic induction ~B, we limit the
required regularity of ~B. This allows the model to include magnetic fields ~B with
strong singularities arising from re-entrant corners, which is not feasible for ~B ∈
H1(Ω)d [4].
We discretize the domain into a shape-regular partition Th of quadrilaterals or
hexahedra {K}. Letting P`(K) be the space of polynomials of degree ` on K and
N`(K) the space of Ne´de´lec vector polynomials of the first kind [14] (with P`−1(K)d ⊂
N`(K) ⊂ P`(K)d), we consider the finite dimensional spaces
V h0 = {~Ch ∈ V0|~Ch|K ∈ N`(K),K ∈ Th}, (2.9a)
V h~q = {~Ch ∈ V~q|~Ch|K ∈ N`(K),K ∈ Th}, (2.9b)
Qh0 = {sh ∈ Q0|sh|K ∈ P`(K),K ∈ Th}. (2.9c)
Then the discrete formulation is as follows: Find ( ~Bh, rh) ∈ V h~q ×Q
h
0 such that
a( ~Bh, ~Ch) + c( ~Bh, ~Ch) + b(~Ch, rh) = 0 (2.10a)
b( ~Bh, sh) = 0 (2.10b)
for all (~Ch, sh) ∈ V h0 ×Q
h
0 .
Let B be the vector containing the coefficients of ~Bh with respect to a basis for
V h0 , and let r be the vector containing the coefficients of rh with respect to a basis for
Qh0 . Then, the finite element solution of the weak formulation (2.10) can be computed
by solving a linear system of saddle point structure,
(
A+N Dt
D 0
)(
B
r
)
=
(
f
0
)
. (2.11)
In this equation, f includes boundary data, A is a discretization of the magnetic
diffusion operator ∇ × ( ηµ∇ × ·), N is a discretization of the magnetic convection
operator ∇× (~u× ·), D is the discrete (negative) divergence operator, and Dt is the
discrete gradient.
We now demonstrate some physical aspects of these equations by applying them
to two deterministic example problems. All simulations throughout the paper are im-
plemented using the deal.II finite element library [1] with first order Ne´de´lec elements
for B and bilinear elements for r (i.e. ` = 1 in (2.9)).
2.1. Example Problem: Hartmann Flow. The Hartmann problem [5] is a
classic two-dimensional test problem modeling the flow of an electrically conducting
fluid through a channel in the presence of an externally applied transverse magnetic
field. We pose this problem on the domain [−0.5, 0.5]2 in the presence of the external
magnetic field ~B = (0, 1). The coupled MHD equations admit the exact analytic
solution
~u =
(
cosh(H/2)−cosh(Hy)
cosh(H/2)−1 , 0
)
(2.12a)
~B =
(
Hν sinh(Hy)−2 sinh(H/2)ycosh(H/2)−1 , 1
)
(2.12b)
for this problem, where H =
√
1
νη is the Hartmann number and ν is the kinematic
viscosity of the fluid. Representative images of the components of ~u and ~B in the
4
−0.5 0 0.50
0.2
0.4
0.6
0.8
1
y
u x
 
 
ν = 10−2ν = 10−1
(a) ux(y)
−0.5 0 0.5−2.5−2
−1.5−1
−0.50
0.51
1.52
2.5
y
B x
 
 
ν = 10−2ν = 10−1
(b) Bx(y)
10−2 10−1
10−8
10−6
10−4
10−2
h
||~ B h−
~ B|| 2
 
 
ν = 10−2ν = 10−1h4
(c) Error
Fig. 2.1: Velocity profiles (a), induced magnetic fields (b), and error || ~Bh − ~B||2 (c)
for the Hartmann problem.
(a) ~u (b) ~Bh with η = 10−3 (c) ~Bh with η = 10−2 (d) ~Bh with η = 10−1
Fig. 2.2: Velocity profile and induced magnetic fields for deterministic MHD eddy
problem.
x-direction are plotted in Figure 2.1 for η = 10−2 and µ = 1 with ν = 10−1 and 10−2.
The plots show that smaller viscosity leads to thinner boundary layers in both ~u and
~B and that the magnitude of the induced magnetic field increases with the viscosity.
This simple example demonstrates that fairly small changes in the velocity field can
lead to large changes in the magnetic field.
We pose a kinematic version of the Hartmann problem by prescribing the velocity
defined by (2.12a) over the domain D and imposing ~B×~n = (0, 1)×~n on the bound-
ary ∂D. The exact solution ~B to this kinematic problem is then given by (2.12b).
Applying the finite element formulation (2.10), we obtain the approximation ~Bh to
~B. To validate the deterministic finite element formulation, we plot the convergence
of the error || ~Bh − ~B||2 as h is refined in Figure 2.1.
2.2. Example Problem: MHD Eddy. The physical effects of the resistivity η
are demonstrated by a two-dimensional benchmark problem considered in [12], which
models the effect of an eddy on a magnetic field. We prescribe the velocity field
~u(x, y) =
(
cos(pix)
pi 32y(1− 4y
2)3
− sin(pix)(1− 4y2)4
)
(2.13)
on the domain Ω = [− 12 ,
1
2 ]
2 and a vertical magnetic field on the boundary with the
condition ~B×~n = (0, 1)×~n on ∂D. Figure 2.2 demonstrates the effect of the resistivity
5
on the induced magnetic field ~B for three values of η on a 64×64 element mesh. In this
figure, we have plotted the velocity profile defined by (2.13) as well as the magnetic
field lines for the solution ~Bh with η = 10−3, 10−2, and 10−1. The figure demonstrates
that two competing physical processes are at play in the kinematics model. First, the
boundary condition corresponds to an external magnetic field applied to the domain.
With infinite resistivity, the velocity field plays no role in the kinematics equations, so
the solution is determined solely by the boundary conditions. In this case, this results
in the uniform vertical magnetic field ~B = (0, 1). For large resistivity, the solution is
dominated by this process. For η = 10−1 for instance, the solution appears to be a
perturbation of the field ~B = (0, 1). The second physical process is governed by the
effect of the velocity field on the magnetic field. For small resistivity, the kinematic
equations are dominated by the convective term ∇× (~u× ~B) which tends to pull the
magnetic field lines in the direction of the velocity field. As the resistivity approaches
zero, the topology of the magnetic field then approaches that of the velocity field. For
this problem, this means that the magnetic field lines should look more like concentric
ellipses as η decreases. This is demonstrated for η = 10−3, where the magnetic field
is nearly “frozen” in the fluid in the center of the domain. Hence, these simulations
show that qualitative characteristics of the magnetic topology can indicate the relative
resistivity of the system. The more the field lines appear to be pulled by the velocity
field (i.e. for this problem, the more swirling in the magnetic field), the smaller the
resistivity.
3. MHD Kinematics with Uncertain Data. Because the equations of MHD
kinematics may involve uncertain quantities, we propose a formulation that incorpo-
rates uncertainty. In particular, we consider the cases where the resistivity η or the
velocity field ~u are uncertain. (Because the permeability µ does not vary significantly
in applications, we will take it as fixed. Equivalently, we can regard any uncertainty in
µ as being absorbed into η.) In this section, we present mathematical representations
for each quantity that incorporate uncertainty by treating η = η(ω) and ~u = ~u(ω) as
random variables. This yields bilinear forms depending on random variables, i.e.
a( ~B, ~C, ω) =
(
η(ω)
µ
∇× ~B,∇× ~C
)
, (3.1a)
c( ~B, ~C, ω) = −
(
~u(ω)× ~B,∇× ~C
)
, (3.1b)
defining a stochastic weak formulation. Thus, by the Doob-Dynkin Lemma, the solu-
tion ( ~Bh, rh) is also a random variable defined on the same sample space [15]. Hence,
each realization of η(ω) and ~u(ω) yields the weak formulation: Find ( ~Bh(ω), rh(ω)) ∈
V h~q ×Q
h
0 such that
a( ~Bh(ω), ~Ch, ω) + c( ~Bh(ω), ~Ch, ω) + b(~Ch, rh(ω)) = 0 (3.2a)
b( ~Bh(ω), sh) = 0 (3.2b)
for all ( ~Ch, sh) ∈ V h0 × Q
h
0 . The solution of any realization of this problem can be
obtained by solving a linear system of the form
(
A(ω) +N(ω) Dt
D 0
)(
B(ω)
r(ω)
)
=
(
f
0
)
. (3.3)
Given this framework, we can employ a Monte-Carlo simulation to obtain statistical
properties of ~Bh. We repeatedly generate independent random instances of η(ω)
6
and ~u(ω) and solve for ~Bh(ω). We then estimate the mean and standard deviation
of ~Bh by the (pointwise) sample mean, which we denote µ( ~Bh)(~x), and the sample
standard deviation, which we denote σ( ~Bh)(~x). A canonical error estimate for Monte-
Carlo simulation [2] states that for N trials, the (pointwise) stochastic error for each
component of µ( ~Bh) satisfies
|E(Bi(~x, ω))− (µ( ~Bh)(~x))i| ≈ 2
(σ( ~Bh)(~x))i√
N
(3.4)
with 95% confidence. Thus, we obtain with 95% confidence the error result
||E( ~B(~x, ω))− µ( ~Bh)(~x)||2 ≈ 2
||σ( ~Bh)(~x)||2√
N
. (3.5)
3.1. Uncertain Velocity. In this section, we consider the case where fluctua-
tions are allowed in the velocity field. We assume that a mean flow, ~u0, is known
and represent the fluctuations by a random variable, ~u∗(ω), with mean zero. We thus
express ~u as the sum of its deterministic and random parts as
~u(ω) = ~u0 + ~u∗(ω). (3.6)
Rather than letting each component of ~u∗ be an independent scalar random variable,
we derive two expressions for ~u∗ from assumptions about physical properties of the
fluid, either that the fluctuation is irrotational or that the fluid is incompressible.
This results in a natural coupling of the components of ~u∗. We explore the effects
of an uncertain velocity field on the MHD kinematics system in Sections 3.1.1 and
3.1.2 by applying fluctuations of these types to the Hartmann flow problem detailed
in Section 2.1.
3.1.1. Test Problem 1: Irrotational Fluctuations. If the random fluctua-
tions of the fluid are irrotational (∇× ~u∗ = 0), then ~u∗ is a conservative vector field
and can be written as the gradient of a scalar potential φ, i.e.
~u∗ = ∇φ. (3.7)
Under this assumption, only the random scalar field φ needs to be specified in order
to define ~u. We assume the potential field to vary continuously and to be spatially
correlated. These assumptions are satisfied if we assume φ to be a stationary random
field with the covariance function defined by
C(~x, ~y) = σ2e−
||~x−~y||2
` . (3.8)
Here, σ2 is the variance and ` is a correlation length. Clearly, the covariance is
greatest when the Euclidean distance between the points ~x and ~y is small. In effect,
this covariance function generates fluctuations in φ on a scale proportional to `.
If φ has mean zero, then φ can be approximated by a truncated Karhunen-Loe`ve
(KL) expansion [7]
φ(~x, ξ) ≈
M∑
i=1
√
λiφi(~x)ξi, (3.9)
7
(a) ux(y) at x = 0 (b) uy(y) at x = 0
Fig. 3.1: Profiles of ux and uy along the line x = 0 for two random instances with
σ2 = 5.0× 10−3, together with the mean profile.
(a) Bx(y) at x = 0 (b) By(y) at x = 0
Fig. 3.2: Profiles of Bx and By along the line x = 0 for two random instances with
σ2 = 5.0× 10−3, together with the deterministic solution obtained from ~u0.
where φi(~x) and λk are the eigenfunctions and eigenvalues of C. We will assume the
random variables {ξi} to be independently and uniformly distributed in [−
√
3,
√
3]M .
We choose M large enough to capture 95% of the total variance [11]; that is,
M∑
i=1
λi > 0.95|D|σ
2. (3.10)
We note that the correlation length affects this requirement, with small ` leading to
large M .
We let the mean velocity profile ~u0 be given by the deterministic Hartmann profile
(2.12a) and introduce fluctuations by letting ~u be defined by (3.6) and (3.9) using the
8
(a) Bx(y) at x = 0 (b) By(y) at x = 0
Fig. 3.3: Profiles of Bx and By for Test Problem 1, plotted along the line x = 0 for
the mean µ( ~Bh) with σ2 = 5.0× 10−3, 6.0× 10−3, and 7.0× 10−3.
covariance function (3.8). We let η ≡ 10−2, ν = 10−2, and ` = 0.1. This correlation
length corresponds to fairly small-scale fluctuations in the velocity field. We compare
the effects of three choices for the variance, σ2 = 5.0×10−3, 6.0×10−3, and 7.0×10−3.
The increase in σ2 corresponds to an increase in the magnitude of the fluctuations.
We present results for this problem discretized on a 64 × 64 element mesh. Sample
random instances of ~u with σ2 = 5.0× 10−3 are plotted in Figure 3.1. In this figure,
the profiles of the components of ~u, ux and uy, are plotted along the line x = 0 and
compared to the mean values (u0)x and (u0)y. The corresponding solutions ~Bh are
plotted in Figure 3.2. These are compared to the deterministic solution obtained with
~u ≡ ~u0. Both the random data ~u and the corresponding solutions demonstrate high
frequency oscillations around the deterministic profiles.
Some results for the Monte-Carlo simulation after 10,000 trials are plotted in
Figure 3.3. In this figure, the profiles of the mean solution µ( ~Bh) are plotted along
the line x = 0 for the three values of σ2. This is compared to the deterministic
solution ~Bh with ~u ≡ ~u0. The Euclidean norm of the pointwise variance of the solution
||σ( ~Bh)(~x)||2 is bounded by 0.28 for σ2 = 5.0×10−3, 0.38 for σ2 = 6.0×10−3, and 0.51
for σ2 = 7.0× 10−3. By (3.5), this implies that the maximum (pointwise) stochastic
error for this problem is approximately 0.01. From Figure 3.3, it can be seen that on
average, irrotational fluctuations in the velocity field result in a slight growth in the
magnitude of the induced magnetic field as compared to the deterministic case. This
growth increases as σ2 increases.
3.1.2. Test Problem 2: Incompressible Flow. If the fluid is assumed to
be incompressible (i.e. ∇ · ~u = 0) and ~u0 is incompressible, then ~u∗ must also be
incompressible (see (3.6). Thus, in this case, we can prescribe ~u∗ to be the curl of a
potential φ. In two dimensions, φ is a scalar, and in three dimensions φ is a vector.
We consider only the 2D case in this study. Then, ~u∗ can be written as
~u∗ =
(
∂φ
∂y
,−
∂φ
∂x
)
, (3.11)
9
(a) Bx(y) at x = 0 (b) By(y) at x = 0
Fig. 3.4: Profiles of Bx and By for Test Problem 2, plotted along the line x = 0 for
the mean µ( ~Bh) with σ2 = 5.0× 10−3, 6.0× 10−3, and 7.0× 10−3.
and the random variable ~u∗ can be computed from the scalar random variable φ. As
above, we let φ be defined by a KL expansion (3.9) with the covariance function C.
Again, we let η ≡ 10−2, ν = 10−2, and ` = 0.1 and consider σ2 = 5.0×10−3, 6.0×10−3,
and 7.0× 10−3.
Mean solution profiles obtained from the Monte-Carlo simulation after 10,000 tri-
als are plotted in Figure 3.4. The normed standard deviation ||σ( ~Bh)(~x)||2 is bounded
by 0.49 for σ2 = 5.0× 10−3, 1.23 for σ2 = 6.0× 10−3, and 3.17 for σ2 = 7.0× 10−3,
corresponding to a maximum stochastic error of about 0.06. Compared to the case of
irrotational fluctuations, both the standard deviation of the solution and the magni-
tude of the induced magnetic field are greater when non-zero vorticity is permitted in
the fluctuations. The difference in the magnitude of the magnetic field is fairly signif-
icant between the two test problems, suggesting that fluid vorticity plays a large role
in generating magnetic fields. Thus, when small-scale rotational behavior is present
in a fluid, simulations based on the mean flow of the fluid may not capture important
magnetic effects.
3.2. Uncertain Resistivity. In this section, we consider the case where just
the resistivity η is a random field over the domain. This is motivated by the fact
that multiple fluids may be present in a physical system and there may be some
epistemic uncertainty in the distribution of the fluids throughout the domain. In
practical applications, the resistivity can range over orders of magnitude between two
fluids. For example, in aluminum electrolysis, the resistivity of liquid aluminum is
approximately 4.0 × 10−3Ωm, while the resistivity of the fluid bath from which the
aluminum is reduced is approximately 2.9× 10−7Ωm [6]. We propose defining
η(~x, ω) = 10β(~x,ω), (3.12)
where β is a random scalar field yet to be specified. This expression both emphasizes
the variability in the order of magnitude of η and guarantees that η > 0. We first
use this to investigate the effects of uncertain resistivity on the MHD eddy problem
10
(a) P1 (b) P2 (c) P3
Fig. 3.5: Domain partitionings considered in Test Problem 3.
(a) η = 10−1 on D1
η = 10−3 on D2
(b) η = 10−3 on D1
η = 10−1 on D2
Fig. 3.6: Instances of ~Bh obtained with partitioning P2 for Test Problem 3.
discussed in Section 2.2, for two different choices of β. We then consider a three-
dimensional extension of the MHD eddy problem
3.2.1. Test Problem 3: Piecewise constant β. In this section, we assume
that the domain is occupied by multiple immiscible fluids with different resistivities.
In this setting, we let β be a piecewise constant scalar field over the domain. We
partition the domain D into n subdomains P = {Dk}nk=1 and let β be a constant on
each of these subdomains. If we assume the resistivity to be uncertain on each of
the subdomains, then we can let β(·, ~ξ)|Dk = ξk where ~ξ = [ξ1, . . . , ξn]
t is a random
vector.
We investigate the effect of a piecewise constant resistivity by considering the
MHD eddy problem, defining ~u by (2.13) on the domain [− 12 ,
1
2 ]
2. We consider three
partitionings of the domain:
P1 = {[− 12 ,
1
2 ]
2}, (3.13)
P2 = {[− 14 ,
1
4 ]
2, [− 12 ,
1
2 ]
2\[− 14 ,
1
4 ]
2}, (3.14)
P3 = {[− 12 , 0)× [0,
1
2 ], [0,
1
2 ]× [0,
1
2 ], [−
1
2 , 0)× [−
1
2 , 0), [0,
1
2 ]× [−
1
2 )}, (3.15)
as shown in Figure 3.5. We let each ξi be independently and uniformly distributed
in the interval [−1.0,−3.0]. Defining η by (3.12), we obtain E(η) ≡ η0 ≈ 2.1 × 10−2
independent of the partitioning. With partitioning P1, the resistivity is a random con-
stant over the domain. Thus, the solutions plotted in Figure 2.2 are representative
instances of magnetic fields that may be induced for this partitioning. With parti-
tioning P2, the resistivity in the center of the domain may differ from the resistivity
11
(a) η = 10−3 on D1, D4
η = 10−1 on D2, D3
(b) η = 10−3 on D2, D3
η = 10−1 on D1, D4
(c) η = 10−3 on D1, D3
η = 10−1 on D2, D4
(d) η = 10−3 on D3, D4
η = 10−1 on D1, D2
Fig. 3.7: Instances of ~Bh obtained with partitioning P3 for Test Problem 3.
(a) ~Bh with η ≡ η0 (b) µ( ~Bh) for P1 (c) µ( ~Bh) for P2 (d) µ( ~Bh) for P3
Fig. 3.8: Deterministic and mean magnetic field lines compared for Test Problem 3.
near the boundaries, and this can result in in different kinds of behavior than seen
for a constant resistivity. Two examples of the kind of behavior we may obtain are
plotted in Figure 3.6. In the first example, the resistivity is much larger in the center
subdomain, and one can see that the behavior of the magnetic field is more charac-
teristic of a larger resistivity in the center. Near the boundary of D, the field lines are
similar to those obtained with η = 10−3 on the entire domain. Around the interface
of the two subdomains, the character of the field lines shifts. The second example
in Figure 3.6 shows the opposite case, where the resistivity is smaller in the center
subdomain. In this case, one can see behavior characteristic of a large resistivity near
the boundary of D and behavior characteristic of a small resistivity in the center.
Figure 3.7 depicts some examples of magnetic fields that can result from partitioning
P3. As with partitioning P2, it can be seen that the character of the field lines in
a particular subdomain are characteristic of the resistivity on that subdomain, and
at subdomain interfaces there is a shift in the character of the magnetic field. Fig-
ures 3.6 and 3.7 demonstrate not only that a discontinuous resistivity can produce
profoundly different solutions than a constant resistivity, but also that the resistivity
in a particular region can be approximated by considering the character of the field
lines in that region.
Some results for the Monte-Carlo simulation after 10,000 trials on a 64×64 element
mesh are shown in Figures 3.8 and 3.9. The field lines for the deterministic case where
η = E(η(ω)) are compared to those obtained from the mean µ( ~Bh)(~x) in Figure 3.8,
and the norm of the standard deviation ||σ( ~Bh)||2 is plotted in Figure 3.9. Note that
12
(a) P1 (b) P2 (c) P3
Fig. 3.9: Euclidean norm of standard deviation ||σ( ~Bh)||2 for Test Problem 3.
the probability that the resistivity is greater than the mean resistivity (Pr(η > η0) ≈
0.33) is less than the probability that the resistivity is less than the mean resistivity
(Pr(η < η0) ≈ 0.67). Despite this, the mean field lines in each example resemble
ones that would arise with η > η0. Because this effect occurs for partitioning P1, it
appears that it is due primarily to the variability of η. When the resistivity is allowed
to vary for this test problem, the mean magnetic field is dominated by the qualities
of magnetic fields induced by larger resistivities. Figure 3.8 shows that the large-
scale behavior of the mean magnetic field is very similar for the three partitionings,
although some small-scale differences are present. For partitioning P2, the behavior
in the center subdomain is consistent with a slightly smaller resistivity than in the
rest of the domain. With partitioning P3, the field lines appear to correspond to
a larger resistivity around the origin where the four subdomains meet. The results
suggest that variability in the resistivity has a more significant effect on the mean
magnetic field than the presence of multiple fluids with different resistivities. Although
multiple resistivities may result in random instances of ~B that differ significantly
from a constant resistivity, as shown in Figures 3.6 and 3.7, the means do not differ
dramatically from solutions for constant resistivities, with the main differences arising
at subdomain interfaces.
Figure 3.9 shows the norm of the standard deviation of ~Bh over the domain.
Because ||σ( ~Bh)(~x)||2 ≤ 0.09, we know from (3.5) that the stochastic error is bounded
as
||E( ~Bh(~x, ω))− µ( ~Bh)(~x)||2 ≤ 2
0.09
√
10, 000
= 1.8× 10−3 (3.16)
with 95% confidence. From the figure, it can be seen that the standard deviation
is affected by the partitioning of the domain. In this case, the standard deviation
increases as the number of subdomains increases. The average standard deviation over
the domain for partitionings P1, P2, and P3 is approximately 2.9 × 10−2, 3.0 × 10−2,
and 4.1 × 10−2. Furthermore, the standard deviation tends to be larger along the
subdomain interfaces.
3.2.2. Test Problem 4: β as a truncated KL expansion. In this section,
we assume that the resistivity of the system varies continuously over space. This
situation may arise when fluids can blend together, such as when different liquid
metals are combined into an alloy. This scenario can be modeled by supplying β with
the spatially correlated covariance function C of (3.8). If β has mean β0(~x), then the
13
(a) ` = 0.1 (b) ` = 1.0 (c) ` = 10.0
Fig. 3.10: Random realizations of η for Test Problem 4.
(a) ~Bh with η ≡ η0 (b) µ( ~Bh) with ` = 0.1 (c) µ( ~Bh) with ` = 1.0 (d) µ( ~Bh) with ` = 10.0
Fig. 3.11: Deterministic and mean magnetic field lines compared for Test Problem 4.
truncated KL expansion for β can be written as
β(~x, ξ) ≈ β0(~x) +
M∑
i=1
√
λiβi(~x)ξi. (3.17)
We assume the random variables {ξi} to be independently and uniformly distributed
in [−
√
3,
√
3]M and choose M large enough to capture 95% of the total variance.
We consider the domain [− 12 ,
1
2 ]
2 and prescribe the velocity field (2.13) with η
defined by (3.12). We let β have mean β0 = −2.0 and variance σ2 = 0.4. This defines
the mean of η to be E(η(ω)) = η0 ≈ 1.5×10−2 (note that η0 6= 10E(β)). We discretize
the problem on a 64×64 element mesh. We consider three choices of correlation length
` = 0.1, 1.0, and 10.0 resulting in truncated KL expansions of length 1834, 39, and
2. Representative realizations of β are plotted in Figure 3.10 for these choices. The
plots demonstrate that as ` increases the resistivity varies over a smaller range for a
given realization. Realizations with ` = 0.1, thus, represent heterogeneous systems in
which fluids with disparate resistivities have not been mixed well. As ` increases the
realizations become more homogeneous in resistivity, corresponding to later stages in
a mixing process. With ` = 10.0, η varies over a very small range in a particular
instance, and in this respect it is similar to the constant resistivity obtained for Test
Problem 3 with partitioning P1. Thus, large correlation lengths correspond to an
uncertain final resistivity at the end of a mixing process.
Figures 3.11 and 3.12 show some results of the Monte-Carlo simulation after
10,000 trials. Figure 3.11 shows the magnetic field lines associated with the mean
14
(a) ` = 0.1 (b) ` = 1.0 (c) ` = 10.0
Fig. 3.12: Euclidean norm of standard deviation ||σ( ~Bh)||2 for Test Problem 4.
of ~Bh with each value of `. These results are compared with the solution to the
deterministic problem for η ≡ η0. It can be seen that with ` = 0.1, the mean field
lines are very similar to the field lines obtained from the deterministic problem. As
` increases, the field lines differ more from the deterministic solution. In fact, the
behavior exhibited for larger ` appears to result from a resistivity greater than the
mean resistivity η0; that is, as ` increases, the magnetic field lines are drawn more
toward the infinite resistivity solution ~B = (0, 1).
Figure 3.12 shows the norm of the standard deviation of ~Bh over the domain.
Because ||σ( ~Bh)(~x)||2 ≤ 0.05, the stochastic error is bounded above by 1.0 × 10−3
with 95% confidence. From this figure, it can be seen that the standard deviation
increases as ` increases. For small `, the mean solution is not only more similar to
the deterministic solution, but the standard deviation is also smaller. This suggests
that the induced magnetic field is more responsive to variations in the resistivity of
homogeneous systems than to small-scale variations in the resistivity of heterogeneous
systems. Although the resistivity can vary over a very large range in a heterogeneous
system, the fluctuations within random instances have little effect on the induced
magnetic field. On the other hand, profound differences from the deterministic case
are seen when the resistivity of a homogeneous system is uncertain. This is consistent
with the results obtained for Test Problem 3, where we found that variability in the
resistivity had stronger effects on the mean solution than the presence of multiple
resistivities.
3.2.3. Test Problem 5: Uncertain Resistivity in 3D. In this section, we
examine the behavior of solutions obtained for three-dimensional models with uncer-
tain piecewise constant resistivity. We take as our domain the cube [− 12 ,
1
2 ]
3, and
consider a generalization of the velocity field (2.13) in which the component in the
z-direction is identically 1,
~u(x, y, z) =


cos(pix)
pi 32y(1− 4y
2)3
− sin(pix)(1− 4y2)4
1

 , (3.18)
which describes a swirling flow in the z-direction. Again, we prescribe a magnetic
field in the y-direction for the boundary condition, ~B × ~n = (0, 1, 0)× ~n.
We define the resistivity by (3.12) using piecewise constant β on the partitionings
P1 = {[− 12 ,
1
2 ]
3}, (3.19)
P2 = {[− 14 ,
1
4 ]
2 × [− 12 ,
1
2 ], ([−
1
2 ,
1
2 ]
2\[− 14 ,
1
4 ]
2)× [− 12 ,
1
2 ]}. (3.20)
15
(a) ~Bh at z = −
1
4 (b)
~Bh at z = 0 (c) ~Bh at z =
1
4
(d) µ( ~Bh) at z = −
1
4 for
P1
(e) µ( ~Bh) at z = 0 for P1(f) µ( ~Bh) at z =
1
4 for P1
(g) µ( ~Bh) at z = −
1
4 for
P2
(h) µ( ~Bh) at z = 0 for P2(i) µ( ~Bh) at z =
1
4 for P2
Fig. 3.13: Deterministic fields lines (top), mean field lines with partitioning P1 (mid-
dle), and mean field lines with partitioning P2 at cross sections z = − 14 (left), z = 0
(middle), and z = 14 (right) for Test Problem 5.
These are extensions of partitionings P1 and P2 from Test Problem 3 in the z-direction.
On each subdomain Dk, we let β(·, ~ξ)|Dk = ξk where ξk is uniformly distributed in
[−1.0,−3.0]. The mean resistivity in this case is η0 ≈ 2.2 × 10−2. Field lines of
the sample mean computed from 1,000 Monte-Carlo simulations on a 16 × 16 × 16
element mesh are compared to deterministic field lines for η ≡ η0 in Figure 3.13.
Field lines are computed in the x-y cross sections at z = − 14 , 0, and
1
4 . From the
deterministic field lines, it can be seen that near the bottom of the domain, the
magnetic field is more dominated by the infinite resistivity solution ~B = (0, 1, 0), but
as z increases the velocity field has a larger effect on the magnetic field. Comparing
to the two-dimensional case, this is as if the resistivity decreases as z increases. We
see the same effect for the mean solution. Unlike for the 2D test problems, the mean
magnetic field looks like a deterministic solution with η < η0. Comparing the mean
field lines obtained for the two different partitionings, differences due to the multiple
resistivities present in partitioning P2 are very slight. As in 2D, this shows that the
16
(a) P1 at z = − 14 (b) P1 at z = 0 (c) P1 at z =
1
4
(d) P2 at z = − 14 (e) P2 at z = 0 (f) P2 at z =
1
4
Fig. 3.14: ||σ( ~Bh)(~x)||2 at cross sections for partitionings P1 (top) and P2 (bottom)
for Test Problem 5.
effects of variability in the resistivity are stronger than the effects of discontinuities
in the resistivity. Cross sections of the norm of the standard deviation are plotted
in Figure 3.14. This shows that variability in the magnetic field is greatest at the
center cross section z = 0. Furthermore, the figure demonstrates that, as in Test
Problem 3, a greater number of subdomains leads to a larger standard deviation that
is distributed throughout more of the domain.
4. Linear Solvers for the Discretized Kinematics System. In Monte-Carlo
simulation, for each realization of the uncertain quantities, a linear system of the form
Ax = b must be solved, where
A = A(ω) =
(
A(ω) +N(ω) Dt
D 0
)
. (4.1)
Because many realizations are required to produce accurate statistical results, it is
imperative that the linear solver be efficient and robust over random variations in
the parameters η and ~u. The linear systems are sparse, nonsymmetric and indefinite,
and, depending on the level of spatial refinement, they can be very large. Hence, a
preconditioned iterative method such as preconditioned GMRES is a natural choice
of solver for these systems. Motivated by the results of [13], many successful block
preconditioners have been developed for solving similar saddle point systems. Follow-
ing this line of research, we develop a generalization of a preconditioner proposed in
[8] for the time-harmonic Maxwell equations to be used for the discretized kinematics
equations.
The system studied in [8] can be considered a special case of (2.5) in which ηµ ≡ 1
and ~u ≡ 0. If we let Aˆ be a discretization of the unscaled magnetic diffusion operator
17
∇×∇×, then the coefficient matrix of the resulting linear system can be written
Aˆ =
(
Aˆ Dt
D 0
)
. (4.2)
The preconditioner developed in [8] is of the form
Pˆ =
(
QB + Aˆ 0
0 Lr
)
, (4.3)
where QB is the mass matrix for B and Lr is a discrete Laplacian on the magnetic
pseudo-pressure space. A generalization of this preconditioner for use with the system
A is
Pk =
(
kQB +A+N 0
0 1kLr
)
, (4.4)
where k > 0 is a constant to be specified.
4.1. Analysis of Eigenvalues. We give a complete analysis of this precondi-
tioner for the case where N ≡ 0. The performance of preconditioner Pk for system A
is governed by the eigenvalues λ of the generalized eigenvalue problem
(
A Dt
D 0
)(
B
r
)
= λ
(
kQB +A 0
0 1kLr
)(
B
r
)
. (4.5)
Defining n = dim(B) and m = dim(r), this has a total of n + m eigenvalues. From
the bottom row of (4.5), we obtain r = kλL
−1
r DB. Substituting this into the top row
of (4.5) gives
λAB + kDtL−1r DB = λ
2(kQB +A)B. (4.6)
Through a discrete Hodge decomposition, B can be written as the sum of its discrete
curl-free part BA and its discrete divergence-free part BD (i.e. B = BA +BD, where
ABA = AtBA = 0 and DBD = 0). Then (4.6) can be rewritten as
λABD + kDtL−1r DBA = λ
2kQB(BA + BD) + λ2ABD. (4.7)
Let the norm induced by a symmetric positive definite matrix X be denoted || · ||X =
〈X·, ·〉1/2. Taking the inner product of (4.7) with BA and using the relations
〈QBBA,BD〉 = 〈QBBD,BA〉 = 0, (4.8a)
〈DtL−1r DBA,BA〉 = ||BA||
2
QB , (4.8b)
proven in [8], we have
k||BA||2QB = λ
2k||BA||2QB . (4.9)
Because there are at least m linearly independent vectors satisfying BA 6= 0, this
means that (4.5) has eigenvalues λ = ±1 each with multiplicity at least m.
Insight into the remaining n−m eigenvalues can be obtained by taking the inner
product of BˆD with (4.7), yielding
λk||BD||2QB = (1− λ)||BD||
2
A. (4.10)
18
From this equation, it is clear that 0 ≤ λ ≤ 1. These eigenvalues can be further
bounded using the discrete coercivity condition. In [8], this condition is written in
terms of the unscaled norm || · ||Aˆ as
||BD||2Aˆ ≥ α
(
||BD||2Aˆ + ||BD||
2
QB
)
, (4.11)
where the constant α ∈ (0, 1) is independent of the mesh. It can be shown that
ηm
µ ||BD||
2
Aˆ
≤ ||BD||2A ≤
ηM
µ ||BD||
2
Aˆ
, (4.12)
where ηm = min~x∈D{η(~x)} and ηM = max~x∈D{η(~x)}, from which we can obtain the
coercivity condition in terms of the scaled norm || · ||A:
||BD||2A ≥ α
(
ηm
ηM
||BˆD||2A +
ηm
µ ||BD||
2
QB
)
. (4.13)
Applying this inequality to (4.10), we obtain
1
1 + kµαηm
ηM−αηm
ηM
≤ λ ≤ 1. (4.14)
The constant ηM−αηmηM is necessarily smaller than 1; thus, we can write the weaker
bound
1
1 + kµαηm
< λ ≤ 1. (4.15)
If we let k = ηmµ , this bound depends only on the coercivity constant α. This depen-
dence on α is similar to that obtained in [8] for Aˆ preconditioned by Pˆ. Because α
is independent of the mesh, letting k = ηmµ defines a preconditioner which should be
robust with respect to both mesh refinement and variations in the resistivity.
When N 6= 0, much of the same analysis applies. Because BA is curl-free, N
satisfies 〈NB,BA〉 = 0. Given this relationship, the presence of N does not affect the
two eigenvalues λ = ±1 with multiplicity m. However, the remaining eigenvalues are
approximated by
(
λ
(
ηm
µ + k
)
− k
)
|BD|2QB = (1− λ)
(
|BD|2A + 〈NB,BD〉
)
. (4.16)
Because 〈NB,BD〉 can become negative, it is difficult to say more about these eigen-
values, but in practice, this preconditioner proves to be effective. This will be demon-
strated in the following section.
4.2. Numerical Results. Because the number of preconditioned GMRES it-
erations required for convergence depends on the linear system obtained from the
random data ~u(ω) and η(ω), we can regard these iteration counts as a random vari-
able. Consequently, we obtain an estimate of the mean of the number of iterations
required for convergence using the sample mean from a Monte-Carlo simulation. The
GMRES iteration continues until the stopping criterion
||b−Ax|| ≤ 10−8||b|| (4.17)
is satisfied. We compute the average number of iterations to reach this criterion over
the Monte-Carlo simulation. We use the preconditioner Pk defined by (4.4) with
19
````````````Test Problem
Mesh
16× 16 32× 32 64× 64
1 (σ2 = 6.0× 10−3) 5.99 5.95 5.83
2 (σ2 = 6.0× 10−3) 5.78 5.77 5.61
3 (P2) 4.58 3.93 3.19
4 (` = 1.0) 4.98 4.32 3.05
4× 4× 4 8× 8× 8 16× 16× 16
5 (P2) 5.16 4.88 4.47
Table 4.1: Preconditioned GMRES iteration counts.
k = ηmµ . Direct methods are used to solve the subsidiary systems corresponding to
the blocks kQB+A+N and 1kLr in Pk. While direct methods are viable for the small
problems investigated in this study, we note that effective multigrid solvers have been
developed for systems similar in form to the block kQB +A+N (see e.g. [9, 10] and
the references therein).
Average iteration counts for each test problem are reported in Table 4.1, each on
three different meshes. These results demonstrate that the preconditioner is highly
effective for both problems with fluctuations in the velocity field (Test Problems 1 and
2) and those with heterogeneous resistivities (Test Problems 3, 4, and 5). Furthermore,
the preconditioner is robust with respect to the mesh, with average iteration counts
decreasing as the mesh is refined.
5. Conclusion. We have presented a numerical method for simulating the kine-
matics of MHD when either the velocity field or the magnetic resistivity of the fluid is
uncertain, applying the method to several steady-state test problems. We have mod-
eled the effect of random perturbations in the velocity field on the induced magnetic
field. In particular, we have demonstrated that, on average, fluctuations with non-zero
vorticity have a large global effect on the induced magnetic field. This supports the
theory that small-scale turbulent flow is necessary for dynamo action. These results
also suggest that simulations based on mean flow may underpredict the magnitude
of magnetic fields. We have also demonstrated that uncertainty in the distribution
of the resistivity can result in different magnetic topologies than in the deterministic
case. Our results show that this effect is most pronounced when the resistivity in
large regions of the domain is uncertain. On the other hand, we have found that,
on average, the induced magnetic field is largely insenstive to small-scale fluctuations
in the resistivity, even when these fluctuations vary over several orders of magnitude
throughout the domain.
In this study, we have introduced several stochastic models for uncertain quanti-
ties in the context of MHD kinematics. By expressing the resistivity as a piecewise
constant field or a truncated KL expansion, we allow for this quantity to be mod-
eled by a vector of independent random variables. We have also proposed physically
motivated expressions for the velocity field that not only couple its components in
a natural way but also require the construction of only one random field in two di-
mensions. We have employed Monte-Carlo simulation to obtain mean and variance
data, but the stochastic expressions introduced here can be used directly in more so-
phisticated uncertainty quantification methods such as stochastic Galerkin, stochastic
20
collocation, and quasi-Monte-Carlo methods.
In addition, we have developed a preconditioner for the discrete kinematics equa-
tions which is robust over variations in both the resistivity and the velocity field. This
is important because many linear systems need to be solved in order to obtain accu-
rate probabilistic distributions of the solution ~Bh. Because this preconditioner is mesh
independent, it allows for the possibility of larger-scale MHD kinematics simulations
in both two and three dimensions. Furthermore, this preconditioner can be useful in
fully coupled MHD models in which the resistivity and velocity field can fluctuate due
to the coupling of the kinematics equations to the Navier-Stokes equations.
REFERENCES
[1] W. Bangerth, T. Heister, G. Kanschat, et al., deal.II Differential Equations Analysis
Library, Technical Reference. http://www.dealii.org.
[2] R. Caflisch, Monte Carlo and quasi-Monte Carlo methods, Acta Numerica, 7 (1998), pp. 1–49.
[3] R. Codina and N. Herna´ndez, Stabilized finite element approximation of the stationary
magneto-hydrodynamics equations, Computational Mechanics, 38 (2006), pp. 344–355.
[4] M. Costabel and M. Dauge, Weighted regularization of Maxwell equations in polyhedral
domains, Numerische Mathematik, 93 (2002), pp. 239–277.
[5] P. Davidson, An Introduction to Magnetohydrodynamics, Cambridge University Press, 2001.
[6] J.-F. Gerbeau, C. Le Bris, and T. Lelie`vre, Mathematical Methods for the Magnetohydro-
dynamics of Liquid Metals, Oxford University Press, 2006.
[7] R. Ghanem and P. Spanos, Stochastic Finite Elements: A Spectral Approach, Dover Publica-
tions, 2003.
[8] C. Greif and D. Scho¨tzau, Preconditioners for the discretized time-harmonic Maxwell equa-
tions in mixed form, Numerical Linear Algebra with Applications, 14 (2007), pp. 281–297.
[9] R. Hiptmair and J. Xu, Nodal auxiliary space preconditioning in H(Curl) and H(Div) spaces,
SIAM Journal on Numerical Analysis, 45 (2007), pp. 2483–2509.
[10] J. Hu, R. Tuminaro, P. Bochev, C. Garasi, and A. Robinson, Toward an h-independent al-
gebraic multigrid method for Maxwell’s equations, SIAM Journal on Scientific Computing,
27 (2006), pp. 1669–1688.
[11] J. Brown Jr., Mean square error in series expansions of random functions, Journal of the
Society of Industrial and Applied Mathematics, 8 (1960), pp. 28–32.
[12] J. Liu, S. Tavener, and H. Chen, ELLAM for resolving the kinematics of two-dimensional
resistive magnetohydrodynamic flows, Journal of Computational Physics, 227 (2007),
pp. 1372–1386.
[13] M. Murphy, G. Golub, and A. Wathen, A note on preconditioning for indefinite linear
systems, SIAM Journal on Scientific Computing, 21 (2000), pp. 1969–1972.
[14] J.-C. Ne´de´lec, Mixed finite elements in R3, Numerische Mathematik, 35 (1980), pp. 315–341.
[15] M. Rao and R. Swift, Probability Theory with Applications, Springer, 2006.
[16] N. Ben Salah, A. Soulaimani, and W. Habashi, A finite element method for magnetohy-
drodynamics, Computational Methods in Applied Mechanics and Engineering, 190 (2001),
pp. 5867–5892.
[17] D. Scho¨tzau, Mixed finite element methods for stationary incompressible magneto-
hydrodynamics, Numerische Mathematik, 96 (2004), pp. 771–800.
21