Vector potential-based MHD solver for non-periodic flows using Fourier continuation
expansions
Mauro Fontana
a
, Pablo D. Mininni
a
, Oscar P. Bruno
b
, Pablo Dmitruk
a
a
Universidad de Buenos Aires, Facultad de Ciencias Exactas y Naturales, Departamento de F ́ısica, & IFIBA, CONICET, Ciudad
Universitaria, Buenos Aires 1428, Argentina.
b
Computing and Mathematical Sciences, Caltech, Pasadena, CA 91125, USA.
Abstract
A high-order method to evolve in time electromagnetic and velocity fields in conducting fluids with non-periodic boundaries
is presented. The method has a small overhead compared with fast FFT-based pseudospectral methods in periodic
domains. It uses the magnetic vector potential formulation for accurately enforcing the null divergence of the magnetic
field, and allowing for different boundary conditions including perfectly conducting walls or vacuum surroundings, two
cases relevant for many astrophysical, geophysical, and industrial flows. A spectral Fourier continuation method is used
to accurately represent all fields and their spatial derivatives, allowing also for efficient solution of Poisson equations
with different boundaries. A study of conducting flows at different Reynolds and Hartmann numbers, and with different
boundary conditions, is presented to study convergence of the method and the accuracy of the solenoidal and boundary
conditions.
Keywords:
MHD; Non-periodic boundary conditions; Fourier continuation; Magnetic vector potential; Direct numerical
simulations.
1. Introduction
Numerical solutions to partial differential equations
(PDEs) have been a cornerstone of engineering and physics
research in the last decades, and their ubiquity has grown
as computational power has increased. Numerical treat-
ment is particularly relevant when dealing with problems
modeled by non-linear PDEs, which are often impossible
to treat analytically even for simple settings. The develop-
ment of algorithms to yield improved numerical solutions,
and to make a more efficient use of computational resources,
remains to the present day an important research topic. In
the particular case of fluid dynamics and plasma physics,
if periodic boundary conditions can be used (e.g., when
considering the bulk flow dynamics), employing Fourier rep-
resentations is a well established method, as it tends to be
optimal when considering the accuracy and computational
efficiency of the numerical scheme [1, 2].
There is, nevertheless, a plethora of physical scenar-
ios that cannot be successfully modeled employing PDEs
with periodic boundary conditions. In those situations
the zoology of numerical techniques is plentiful, and the
method of choice is usually problem dependent [
3
,
4
,
5
,
6
,
7
,
8
,
9
,
10
,
11
]. One strategy among many possible is to
employ Fourier representations and adjust boundary condi-
tions, most commonly, using either penalization methods
[
12
,
13
], or periodic extensions of the non-periodic fields
Email address:
mfontana@df.uba.ar
(Mauro Fontana)
[
14
,
15
]. However, the former tend to result in low or-
der methods, whereas the latter can become prohibitively
computationally expensive and hence, historically, Fourier-
based methods have not seen much use for non-periodic
problems, resulting in the use of other bases (in many cases
more expensive to transform), or of lower-order methods
for many applications.
In recent years, an efficient Fourier representation for
fields with arbitrary boundary conditions was introduced,
which uses Gram polynomials (sometimes referred to as
discrete Chebyshev polynomials) to obtain high-order peri-
odic extensions at a marginal computational cost, resulting
in the so-called FC-Gram method [
16
,
17
]. Besides yield-
ing high order accuracy, this technique produces spatially
dispersionless derivatives [
18
] — in view of its reliance on
Fourier expansions — and hence phase speed propagation
errors might arise solely as a result of the time stepping
strategy. Another advantage of the FC-Gram method is
that Poisson equations, which are common in some for-
mulations of the incompressible hydrodynamic equations,
can be easily and efficiently solved in bounded domains
[
19
]. As a result, to the present day the FC-Gram method
has been successfully employed in a wide range of PDE
problems [17, 20, 21, 19].
Magnetohydrodynamic (MHD) equations, and related
equations of electromagnetism and plasma physics, have
many features in common with hydrodynamics but also
their own complexities, and as a result have remained
so far unexplored using FC-Gram methods. These non-
linear equations, which model conducting fluids and the low
Preprint submitted to Elsevier
February 3, 2022
arXiv:2107.07077v2 [physics.comp-ph] 2 Feb 2022
frequency regime of plasmas [
22
], are of major relevance to
numerous research topics such as planetary and celestial
magnetic fields [
23
,
24
,
25
,
26
,
27
,
28
], stellar atmospheres
[
29
,
30
], non-mechanical pumps [
31
,
32
], and fusion research,
to name a few. However, boundary conditions for the
electromagnetic fields can be numerically harder to deal
with, and the magnetic field must be kept solenoidal or
otherwise measurable differences in the numerical solutions
can be found [
33
,
34
]. Moreover, the MHD equations
present a purely magnetic invariant, the magnetic helicity,
whose definition involves the magnetic vector potential [
35
].
It is, therefore, also important for MHD research to be able
to confidently rely on numerical methods which produce
an accurate vector potential and, by proxy, are reliable
for computing all the known ideal invariants of the system
including the magnetic helicity. As a result, methods which
evolve equations for the vector potential instead of the
magnetic field are of particular relevance, as they have the
benefits of explicitly enforcing the solenoidal condition on
the magnetic field and of allowing direct computation of
the magnetic helicity.
In this paper we present an algorithm (and a publicly
available numerical code) for the solution of the incompress-
ible MHD equations in the vector potential formulation,
for a single non-periodic dimension along which FC-Gram
representations are employed to evaluate derivatives. In
particular, the highly relevant scenarios where the bound-
ary conditions correspond to either perfectly conducting
walls or vacuum surroundings are considered, as often
used in simulations of the geodynamo or in industrial ap-
plications. Use of an explicit time integration scheme is
proposed, as we will focus on the moderate and low diffu-
sivity regime, where turbulent behaviour takes place and
hence the Courant–Friedrichs–Lewy (CFL) constraint is
dominated by field advection (i.e., by non-linearities). As
indicated above, for simplicity, the method presented in
this work is restricted to cases with a single non-periodic
direction. It is however straightforward to extend the
method presented here to additional non-periodic dimen-
sions and to general non-rectangular domains, as well as to
the compressible case and to other physical systems such
as Hall-MHD flows, two-fluid plasma descriptions, or other
MHD problems including systems with thermal convection
or in rotating frames. Convergence of the method and
the accuracy in satisfying the solenoidal and boundary
conditions will be illustrated considering the Hartmann
flow.
The main results are: (1) The method retains all the
advantages of spectral methods for dealing with strongly
non-linear systems, including being dispersionless and fast
(as its based on FFTs). (2) Being a high-order method,
and based on a high-order representation of the vector
potential, it satisfies the solenoidal condition on the mag-
netic field down to an error of the order of floating point
arithmetic limitations. (3) Solving Poisson equations (e.g.,
to impose the Coulomb gauge) is easy and well conditioned.
(4) The proposed method can accommodate different types
of boundary conditions for the magnetic field, and in par-
ticular, cases relevant for geophysical, astrophysical and
industrial flows are explicitly treated. And finally, (5) a
modified FC-Gram method allows for flexible and efficient
numerical treatment of boundary conditions on second
derivatives of the fields, and on boundary conditions pre-
scribed by differential equations as in the case of Robin
boundary conditions.
The paper is organized as follows. In Section 2 we intro-
duce the equations and geometry of interest. In Section 3
we derive appropriate boundary conditions for the vector
potential. In Section 4, an overview of the standard FC-
Gram method is presented, and appropriate generalizations
for the boundary conditions in question are discussed. In
particular, we point the reader to Section 4.3 which intro-
duces the modified FC-Gram approach for Robin boundary
conditions — without which the applicability of the FC-
Gram method in this context would be compromised, in
view of its ubiquity in electromagnetic problems and of its
potentially high computing costs. Then, in Section 5 a full
time stepping algorithm is proposed for evolving the MHD
equations. As an application, in Section 6 we consider
the Hartmann flow scenario from a physical standpoint,
while in Section 7 the numerical performance of the algo-
rithm is evaluated. Finally, the conclusions are presented
in Section 8.
2. Governing equations
The standard field formulation for the MHD approxima-
tion of an incompressible conducting fluid is encompassed
in the set of equations
∇
·
u
= 0
,
(1)
∂
u
∂t
+ (
u
·
∇
)
u
=
−
∇
p
+
j
×
b
+
ν
∇
2
u
,
(2)
∇
·
b
= 0
,
(3)
∂
b
∂t
=
∇
×
(
u
×
b
) +
η
∇
2
b.
(4)
Here
u
,p,
b
, and
j
=
∇
×
b
are the velocity, the pressure
per unit mass density, the magnetic field, and the current
density, respectively. The kinematic and magnetic diffusiv-
ities are
ν
and
η
, and the mass density is assumed to be
homogeneous. The magnetic field
b
is written in velocity
(Alfv ́enic) units. The magnetic permeability
μ
is equal to
1, a reasonable choice in most optical, geophysical, and
astrophysical contexts. In this approximation the electric
field
E
is reduced to a secondary role as a result of the
quasineutrality hypothesis, although it can be recovered a
posteriori from the non-relativistic Ohm’s law
η
j
=
E
+
u
×
b
.
(5)
Equation (3) implies that the magnetic field is always
solenoidal, and imposes a condition on advection and
stretching of magnetic field lines [
22
]. In low-order meth-
ods, or in approximate methods, it has been shown that
2
solutions are sensitive to the way this condition is enforced,
or on how well it is satisfied [
33
]. In particular, measurable
differences in physical quantities were reported for different
methods [
34
]. When treating problems analytically, the
condition
∇
·
b
= 0 can be enforced by representing the
magnetic field via the vector potential
a
which is defined
by the relation
∇
×
a
=
b
. This latter expression, how-
ever, does not define a unique vector potential, as adding
a gradient to
a
results in the same magnetic field. To
properly define a vector potential both
∇
·
a
and boundary
conditions compatible with those for the magnetic field
must be prescribed. In non-relativistic MHD, the most
common choice is to use the Coulomb gauge,
∇
·
a
= 0.
Besides guaranteeing
∇
·
b
= 0, the vector potential is
also important to define the magnetic helicity density,
h
m
=
a
·
b
. The volumetric integral of this quantity measures
the topological complexity of magnetic field lines, and
is the only purely magnetic ideal invariant of the MHD
equations, the other two known invariants involving also
the velocity field being the volumetric integral of the total
energy density,
e
=
u
2
+
b
2
, and the volumetric integral of
the cross helicity density,
h
c
=
u
·
b
[
22
,
36
]. It is therefore
desirable in the numerical study of MHD to ensure an
accurate computation of the magnetic helicity, for which
the vector potential is needed.
If the magnetic field is represented via a vector potential
in the Coulomb gauge, integrating Eq. (4) results in the
following set of equations for the MHD approximation of
an incompressible flow in terms of
a
,
∇
·
u
= 0
,
(6)
∂
u
∂t
+ (
u
·
∇
)
u
=
−∇
p
−∇
2
a
×
(
∇
×
a
) +
ν
∇
2
u
,
(7)
∇
·
a
= 0
(8)
∂
a
∂t
=
u
×
(
∇
×
a
) +
η
∇
2
a
−
∇
φ.
(9)
where we used
j
=
−∇
2
a
. The scalar potential
φ
results
from the integration of Eq. (4), and must fulfill the gauge
condition
∇
·
a
= 0. A close look at Eq. (5) reveals that
φ
is no other than the standard electric scalar potential,
recovering the usual expression for the electric field
E
=
−
∇
φ
−
∂
a
∂t
.
(10)
It should also be noted that, under this standard formula-
tion of the incompressible MHD equations, both the pres-
sure and the electric potential must fulfill seldom Poisson
equations, since taking the divergence of Eqs. (7) and (9)
and considering the solenoidality condition yields
∇
2
p
=
−
∇
·
[
(
u
·
∇
)
u
]
,
(11)
∇
2
φ
=
∇
·
[
u
×
(
∇
×
a
)
]
.
(12)
The numerical method presented in this work can com-
pute solutions for Eqs. (7), (9), (11) and (12) after bound-
ary conditions for the physical fields are provided. We will
consider a (0
,
0
,
0)
×
(
L
x
,L
y
,L
z
) cuboid domain, and a
uniform spatial discretization in each direction consisting
of
N
x
×
N
y
×
N
z
gridpoints (see Fig. 1).
3. Boundary conditions for the vector and scalar
potentials
3.1. Perfectly conducting boundary conditions
A relevant scenario for astrophysical and geophysical
applications is that of a magnetofluid that is confined
by a perfectly conducting medium. This is the situation
encountered, e.g., in many simulations of the geodynamo
when treating the boundary with the inner core [
25
]. In
this case the electric field must vanish inside the conductor,
and there cannot exist a time-varying magnetic field in its
interior. Considering now a stationary conductor (
u
=
0
),
the continuity of the normal component of the magnetic
field, and Ohm’s law, the magnetic boundary conditions
that the system must fulfill are
j
|
∂
Ω
×
ˆ
n
=
0
,
(13)
∂
t
b
|
∂
Ω
·
ˆ
n
=
0
,
(14)
where
∂
Ω denotes the boundary surface and
ˆ
n
is a normal
unit vector at
∂
Ω. Note that if
b
·
ˆ
n
= 0 at
t
= 0, then
b
·
ˆ
n
= 0
∀
t >
0. For simplicity we consider this latter
scenario from here onwards.
Assuming now periodic boundary conditions in the
ˆ
x
and
ˆ
y
directions and the presence of a conducting wall
at
z
= 0, it follows that a vector potential satisfying the
Coulomb gauge constraint,
a
|
z
=0
×
ˆ
z
= 0, and (
∂
2
a
/∂z
2
)
|
z
=0
×
ˆ
z
= 0 at
t
= 0 verifies the boundary conditions given by
Eqs. (13) and (14). Note that the normal component of the
vector potential,
a
⊥
=
a
z
, is not directly involved in the
determination of any of the physical boundary conditions.
For this to be the case, however, the Coulomb gauge must
be maintained at all times and thus it is reasonable for
a
z
at the boundary to explicitly enforce the interface gauge
condition, resulting in the set of conditions
a
|
z
=0
×
z
= 0
,
(15)
∂
2
a
∂z
2
∣
∣
∣
∣
∣
z
=0
×
ˆ
z
=
0
,
(16)
∂
a
∂z
∣
∣
∣
∣
z
=0
·
ˆ
z
= 0
,
(17)
φ
|
z
=0
= 0
,
(18)
which guarantee the physical boundary conditions are main-
tained at all times. Even more, choosing
φ
= 0 in the
boundary is in full agreement with the fact that the per-
fect conductor must be an equipotential surface. Although
the boundary conditions above appear at first sight to
overdetermine the problem, it should be noted that some
choices are made using the gauge freedom, and that if
they are enforced exactly the dynamical equations imply
3
Figure 1:
a)
:Illustration of the explored geometry, a cuboid domain, discretized employing a uniform grid in each of
ˆ
x
,
ˆ
y
, and
ˆ
z
directions. For
simplicity, the only non-periodic direction is
ˆ
z
. The actual physical domain is represented in light blue, whereas the spaced used to compute
the FC-Gram continuations is shaded in red.
b)
Sketch of the magnetic field lines inside the domain when considering perfectly conducting
boundary conditions, together with the appropriate boundary conditions.
c)
Sketch of the magnetic field inside and outside the domain
occupied by the fluid for the case of vacuum boundary conditions.
that
a
‖
|
z
=0
=
0
acts only as an initial condition (with the
subindex indicating the components parallel to the bound-
ary), and that together with (
∂
2
a
‖
/∂z
2
)
|
z
=0
=
0
, which
is the physical boundary condition, suffice to enforce the
conditions on the magnetic field and the current density.
3.2. Vacuum boundary conditions
Another relevant astrophysical and geophysical scenario
is that of a magnetofluid which is surrounded by vacuum.
As mentioned before, we assume for this work that in the
magnetofluid
μ
= 1, and hence the magnetic field at the
interface must be continuous, that is
[
b
]
∂
Ω
=
0
=
⇒
b
I
|
∂
Ω
=
b
II
|
∂
Ω
,
(19)
where
I
and
II
denote quantities in the fluid and vacuum
media respectively,
∂
Ω is the boundary between them, and
[
q
]
S
denotes the jump in field
q
across
S
. Expanding
both fields in terms of the respective vector potentials
b
i
=
∇
×
a
i
it is straightforward to get the corresponding
bulk equations for
a
II
in the vacuum domain, as the current
density must necessarily vanish there. Due to the absence
of electric charge, the vacuum electric potential
φ
II
must
be harmonic. The respective equations are
∇
2
φ
II
= 0
,
(20)
∇
2
a
II
= 0
,
(21)
∇
·
a
II
= 0
,
(22)
where, as before, the Coulomb gauge was chosen for the
vector potential.
The physical boundary condition, Eq. (19), implies that
a
must obey
(
∇
×
a
I
)
∂
Ω
= (
∇
×
a
II
)
∂
Ω
,
(23)
whereas regularity conditions on
b
require
a
to be continu-
ous at the boundary, that is,
a
I
=
a
II
.
To address the physical boundary conditions for the
scalar potential, we start by formulating appropriate bound-
ary conditions for the electric field, which are
E
I
×
ˆ
n
|
∂
Ω
=
E
II
×
ˆ
n
|
∂
Ω
,
(24)
E
I
·
ˆ
n
|
∂
Ω
=
I
II
E
II
·
ˆ
n
|
∂
Ω
,
(25)
with
i
the electrical permittivity of medium
i
, and
ˆ
n
a nor-
mal unit vector pointing from
I
towards
II
. For a sizeable
class of conducting fluids the permittivity can be approxi-
mated by that of vacuum, and hence the normal component
of the electric field is continuous in the boundary, which
expressed in terms of the potentials reads
(
∂
a
I
∂t
·
ˆ
n
+
∂φ
I
∂n
)
∂
Ω
=
(
∂
a
II
∂t
·
ˆ
n
+
∂φ
II
∂n
)
∂
Ω
.
(26)
Considering the continuity of the vector potential across
the boundary, the latter equation reduces simply to the
continuity of the normal derivative of the scalar potential.
Analyzing a similar equation for the parallel electric field
leads to the continuity of the electric potential itself.
For the case of periodic boundary conditions in
ˆ
x
and
ˆ
y
, assuming vacuum is present in the semispace
z >
0
and considering only bounded solutions for
z
→ ∞
, the
harmonic equations for
a
II
and
φ
II
admit the solutions
a
II
(
x,y,z
) =
∞
∑
m
=
−∞
∞
∑
n
=
−∞
ˆ
a
II
nm
e
i
(
k
x
n
x
+
k
y
m
y
)
−
γ
nm
z
,
(27)
φ
II
(
x,y,z
) =
∞
∑
m
=
−∞
∞
∑
n
=
−∞
ˆ
φ
II
nm
e
i
(
k
x
n
x
+
k
y
m
y
)
−
γ
nm
z
,
(28)
where
k
x
n
(resp.
k
y
m
) is the
n
-th (resp.
m
-th) wavenumber
in the
ˆ
x
(resp.
ˆ
y
) direction and
γ
nm
= [(
k
x
n
)
2
+ (
k
y
m
)
2
]
1
/
2
.
For the vector potential, the condition
a
z
II
nm
=
i
γ
mn
(
k
x
n
a
x
II
nm
+
k
y
m
a
y
II
nm
)
,
(29)
corresponds to enforcement of the Coulomb gauge.
On the basis of this solution boundary conditions for
the vector potential
a
I
can be obtained easily on the base of
its continuity, the physical boundary conditions in Eq. (23)
and the gauge continuity across the boundary, leading to
the set of conditions
(
d
a
I
nm
d
z
+
γ
nm
a
I
)
∣
∣
∣
∣
∣
∣
z
=0
= 0
,
(30)
4
which represents an homogeneous Robin boundary condi-
tion for each wavenumber component of the vector potential
in the fluid medium. In a similar way, the continuity of
both
∂
z
φ
and
φ
itself leads to a fully analogous Robin
condition
(
d
φ
I
nm
d
z
+
γ
nm
φ
I
)
∣
∣
∣
∣
∣
∣
z
=0
= 0
.
(31)
It is worth emphasizing that these boundary conditions
assume vacuum is present in the semispace
z >
0. Nat-
urally, if vacuum is instead in the semispace
z <
0, the
corresponding harmonic solution has a factor
e
γ
nm
z
(note
the absence of the minus sign), and the resulting boundary
conditions for this case would be
(
d
d
z
−
γ
nm
)
∣
∣
∣
∣
∣
z
=0
= 0
,
(32)
for
φ
and each cartesian component of
a
.
Robin boundary conditions are common in many elec-
tromagnetic problems, and as a result, we now present a
generalization of the FC-Gram method to deal with such
conditions, as well as with the boundary conditions dis-
cussed before for the case of a perfect conductor.
4. Generalization of FC-Gram to boundary condi-
tions in MHD
4.1. FC-Gram with Dirichlet boundary conditions
In order to have an accurate but computationally ef-
ficient representation of the fields and their derivatives,
all the relevant variables are projected onto a Fourier rep-
resentation basis. As both the
ˆ
x
and
ˆ
y
directions have
periodic boundary conditions, the transformation to the
wavenumber domain can be directly computed via standard
FFT operations. That is not the case, however, for the
non-periodic
z
dimension, as the well known Gibbs ringing
phenomenon [
37
] would severely degrade accuracy in the
representation of the fields’ derivatives. To circumvent this
limitation we employ a computationally efficient continu-
ation methodology, known as FC-Gram, first introduced
in [
38
,
16
] and recently utilized in [
19
] for a hydrodynamic
Navier-Stokes solver.
The idea behind a Fourier Continuation (FC) technique
can be summarized as follows for the one dimensional case.
Let
f
(
z
) be a non-periodic function defined over the dis-
crete grid
z
i
=
i
∆
z
, with ∆
z
the (uniform) grid spacing
and
i
= 0
, ... N
−
1, resulting in
f
=
[
f
(
z
0
)
,...,f
(
z
N
−
1
)
]
.
The method generates an efficient continuation spanning
the domain
z
N
, ... ,z
N
+
C
−
1
, where
C
is the number
of points used for the extension operation and is a pa-
rameter of the method. The resulting quantities
f
c
=
[
f
(
z
N
)
, ... ,f
(
z
N
+
C
−
1
)
]
can then be appended to the
values of the function over the original grid, as depicted
in Fig. 2, and a Fourier representation can be computed
for
f
∪
f
c
over the interval [
z
0
,z
N
+
C
). Any derivatives
estimated in the wavenumber domain can then be inverse
Fourier transformed to get accurate representations over
the original [
z
0
,z
N
−
1
] grid.
To obtain
f
c
in a computationally efficient manner, the
FC-Gram method uses only information near the boundary
to estimate
f
c
. To this end, lets first consider the simpler
case where we have
d
values at the end of the domain
f
e
dir
=
[
f
(
z
N
−
d
)
, ... ,f
(
z
N
−
1
)
]
, and we wish to compute
a continuation that smoothly transitions from
f
(
z
N
−
1
) to 0.
The subindex “dir” is used to denote a Dirichlet boundary
condition on the endpoint, i.e., that
f
(
z
N
−
1
) is known.
The resulting extension will hence have a
d
-th order of
approximation at the boundary. Note
d
(the number of
points near the boundary used to compute the continuation)
is an additional parameter of the FC-Gram technique.
Fixing
d
defines a set of discrete polynomials orthogonal
with respect to the inner product
〈
g
|
h
〉
=
d
−
1
∑
i
=0
g
(
x
i
)
h
(
x
i
)
,
(33)
i.e., the Gram polynomials. Note that
x
i
is used in this
equation as the inner product definition is unrelated to
the grid
z
i
introduced before. An operator of length
C
that smoothly transitions each Gram polynomial to zero
(i.e., a blend to zero operator) can be numerically obtained,
resulting in the
C
×
d
blend-to-zero matrix
A
dir
[
20
]. The
subindex “
dir
” again denotes the fact that this matrix is
for Dirichlet boundary conditions. Similarly, an operator
that projects arbitrary function values at
d
grid points onto
each of the Gram polynomials can be obtained, resulting
in a certain matrix
Q
dir
(see [20, 19] for more details).
The computation of operators
A
dir
and
Q
dir
must be
performed utilizing arbitrary precision linear algebra rou-
tines, as it requires decomposing ill-conditioned matrices.
However, it is worth noting that these computations, which
require only a few minutes in a single modern CPU core
(although its computation can be parallelized), must be
performed only one time for a given selection of
C
and
d
, with the resulting operators being utilized henceforth.
Additional details about obtaining these operators can be
found in [
38
,
16
], and also in [
19
] for a case with a very
similar scope to the one proposed in this work.
With the
A
dir
and
Q
dir
operators at hand, finding a
continuation that smoothly transitions
f
to 0 is straight-
forward, requiring only a projection of the boundary end
values
f
e
dir
onto the Gram space, and blending each of those
polynomials to zero. The resulting smooth extension to
zero
̃
f
c
is computed as
̃
f
c
=
A
dir
Q
dir
f
e
dir
.
(34)
It is now possible to consider the original problem in
which we want a smooth transition from
f
e
dir
to the
d
values
at the beginning of the domain
f
b
dir
=
[
f
(
z
0
)
,...,f
(
z
d
−
1
)
]
,
to get a periodic extension. This can be easily solved now by
adding to
̃
f
c
a set of values which smoothly transition from
5
0 to
f
(
z
0
). A suitable continuation which transitions from
f
(
z
N
−
1
) to
f
(
z
0
) and which has
d
−
1 smooth derivatives
at the endpoints can thus be constructed as
f
c
=
A
dir
Q
dir
f
e
dir
+
A
‡
dir
Q
Π
dir
f
b
dir
,
(35)
where
‡
and
Π
denote the row-reversing and column-reversing
operations, respectively. The underlying concept in Eq. (35),
which can be easily appreciated in Fig. 2, is the superposi-
tion of values which smoothly transition from
f
(
z
N
−
1
) to
zero (given by
A
dir
Q
dir
f
e
dir
), and from zero to
f
(
z
0
) (pre-
scribed by
A
‡
dir
Q
Π
dir
f
b
dir
), resulting in the overall smooth
periodic extension.
It is worth pointing out that, as mentioned earlier,
A
dir
and
Q
dir
depend only on the choice of
d
and
C
, and not
on the data to be continued
f
b
dir
,
f
e
dir
. This means that
the ill-conditioning usually found when trying to compute
suitable periodic extensions is encapsulated in the estima-
tion of those matrices. This is specially important for 3D
problems, as the fastest way to obtain the values on the
extended volume (as depicted in Fig. 1) is to simply apply
Eq. (35) to each
xy
-line (i.e., by computing
N
x
×
N
y
1D
continuations). Another detail worth noting is that the
same
A
dir
and
Q
dir
operators can accurately continue both
the real and imaginary parts in case the vectors
f
b
dir
,
f
e
dir
are
complex valued, allowing to commute FFTs and FC-Gram
operations (which is important for efficient parallelization
of the method, see [
39
,
19
]). It is also worth mentioning
that, depending on the problem at hand, great accuracy is
obtained when considering
d
in the range of 4 to 10 and
C
in the range of 15 to 33. This has the added benefit that
consequently both
A
dir
and
Q
dir
fit easily in the L1 cache
of modern CPU cores.
Within the context of PDE solvers, the FC-Gram pic-
ture described above is useful when the function values are
known everywhere, that is, both at every interior point as
well as at the boundary. For the case of having the normal
derivative prescribed at one end of the domain, a solution
within the framework of FC-Gram which is both accurate
and efficient was introduced in [
20
]. A treatment for the
case in which second normal derivatives are prescribed at
the boundary is introduced next in Section 4.2. A straight-
forward extension of these methods to Robin boundary
conditions is not efficient, as in electromagnetism a large
number of such conditions may arise simultaneously with
different Robin coupling parameters for each mode (as, e.g.,
in Eq. 30). Thus, a new FC-Gram algorithm for the Robin
problem that allows quick computation for different Robin
coupling parameters is presented in Section 4.3.
4.2.
Modified FC-Gram method for prescribing the second
normal derivative
We are once again interested in finding a suitable blend
to zero procedure using the information at the end of the
domain, but now for the case where
d
−
1 values at the
interior points and the second derivative at the last grid
point are known. Once found, that blend to zero procedure
0
.
0
0
.
2
0
.
4
0
.
6
0
.
8
1
.
0
−
4
−
3
−
2
−
1
0
1
2
a)
0
-2.5×10
6
0.0
2.5×10
6
1
0.93355
0.93357
0.93359
f
f
c
rob
f
c
neu
2
f
b
f
e
AQ
f
e
A
‡
Q
Π
f
b
0
.
0
0
.
2
0
.
4
0
.
6
0
.
8
1
.
0
−
50
0
50
b)
f
′
f
′
rob
f
′
neu
2
0
.
0
0
.
2
0
.
4
0
.
6
0
.
8
1
.
0
z
10
−
5
10
−
3
|
f
′
−
f
′
|
c)
f
′
rob
f
′
neu
2
4th order FD
Figure 2: Example of a 4
th
order (
d
= 5) FC-Gram procedure applied
to the function
f
(
z
) =
J
0
(25
z
)
e
3
z
−
1 for
N
= 113 grid points in the
[0
,
1] domain and
C
= 15 continuation points. Panel
a)
shows the
function over the discrete grid (in blue) together with the periodic
continuations computed using the second normal derivative
f
c
neu
2
(in green with
×
markers) and the Robin condition
f
′
+ 5
f
,
f
c
rob
(in orange with + markers), as the boundary condition at both
endpoints. Also displayed are the matching values used to perform
the continuations,
f
b
(in yellow circles) and
f
e
(in pink squares), as
well as the blend to zero continuations associated to
f
e
(
AQ
f
e
) with
a red dotted line, and
f
b
(
A
‡
Q
Π
f
b
) with a red dashed-dotted line.
Insets show
f
(0) and
f
(1) together with the values reconstructed
from the boundary conditions. Panel
b)
shows the exact derivative
f
′
(in blue using circle markers) together with the ones computed
spectrally from
f
c
neu
2
(in green with
×
markers) and
f
c
rob
(in orange
with + markers) over the [0
,
1] grid. Panel
c)
shows the absolute
error in the derivative estimation for both
f
c
neu
2
(in green with
×
markers) and
f
c
rob
(in orange with + markers). For comparison, a
same order finite difference error is also exhibited (with a dashed
black line). Note that for display purposes marker spacing is larger
than grid spacing in the whole figure, except for the quantities
f
b
and
f
e
for which grid and marker spacing are the same.
6
can be easily employed to compute a periodic extension,
as it was previously discussed. Therefore, in a similar way
to the ideas introduced for the Dirichlet case, one possible
strategy to deal with this problem within the framework of
the FC-Gram method is to project the last
d
−
1 interior
values plus the prescribed second normal derivative onto
a space of modified Gram polynomials, ones which are
orthogonal with respect to the inner product
〈
g
|
h
〉
=
g
′′
(
x
d
−
1
)
h
′′
(
x
d
−
1
) +
d
−
2
∑
i
=0
g
(
x
i
)
h
(
x
i
)
.
(36)
Analogously to the procedure described in [
16
,
20
], given
an arbitrary uniform grid
x
0
,...,x
d
−
1
with spacing ∆
x
, a
set of such modified Gram polynomials can be obtained
employing the modified Vandermonde matrix
P
neu
2
=
1
x
0
x
2
0
...
x
d
−
1
0
1
x
1
x
2
1
...
x
d
−
1
1
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1
x
d
−
2
x
2
d
−
2
...
x
d
−
1
d
−
2
0
0
2
...
(
d
−
1)(
d
−
2)
x
d
−
3
d
−
1
,
(37)
via its
QR
decomposition
P
neu
2
=
Q
neu
2
R
neu
2
(the subindex
“neu
2
” here stands for the case with Neumann conditions on
the second derivative). It is then possible to proceed and
construct blend-to-zero operators for this set of polynomi-
als in an analogous way to the Dirichlet case, obtaining an
operator
A
neu
2
. From these matrices a blend-to-zero contin-
uation can be computed given interior values and the respec-
tive second derivative,
n
=
[
f
(
x
0
)
, ... ,f
(
x
d
−
2
)
,f
′′
(
x
d
−
1
)
]
,
directly as
A
neu
2
Q
neu
2
n
.
One comparably accurate and more computationally
convenient procedure is to obtain the function value at
x
d
−
1
employing the known function values at
x
0
,...,x
d
−
2
and the second derivative at the endpoint. Once
f
(
x
d
−
1
)
is known, continuation values can be computed using the
Dirichlet blend to zero operator
A
by means of Eq. (35).
For this purpose it is useful to first introduce the standard
Vandermonde matrix
P
dir
=
1
x
0
x
2
0
... x
d
−
1
0
1
x
1
x
2
1
... x
d
−
1
1
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1
x
d
−
2
x
2
d
−
2
... x
d
−
1
d
−
2
1
x
d
−
1
x
2
d
−
1
... x
d
−
1
d
−
1
,
(38)
with an associated QR factorization
P
dir
=
Q
dir
R
dir
. The
corresponding function end value can be obtained from the
relations
P
neu
2
c
=
n
T
,
(39)
P
c
=
d
T
,
(40)
with
c
the coefficient of each monomial term on the respec-
tive Gram basis and
d
=
[
f
(
x
0
)
, ... ,f
(
x
d
−
2
)
,f
(
x
d
−
1
)
]
.
Solving for
c
using the
QR
factorizations of the respective
P
matrices results in
Q
dir
R
dir
R
−
1
neu
2
Q
T
neu
2
n
T
=
d
T
,
(41)
and hence the end value
f
(
x
d
−
1
) is given by
f
(
x
d
−
1
) =
̃
q
neu
2
n
T
,
(42)
where
̃
q
neu
2
is the last row of
Q
dir
R
dir
R
−
1
neu
2
Q
T
neu
2
, that is
( ̃
q
neu
2
)
i
= (
Q
dir
R
dir
R
−
1
neu
2
Q
T
neu
2
)
d
−
1
,i
.
It should be noted that the order of accuracy in the
reconstruction of
f
(
x
d
−
1
) is still
d
, independently of the fact
that a second derivative is prescribed instead of the value of
f
at
x
d
. Thus, if the reconstruction is performed utilizing
d
neu
2
values at the end of the domain (
d
neu
2
−
1 function
values and the prescribed second derivative) after which
a
C
×
d
dir
Dirichlet blend to zero operator
A
dir
is used to
compute continuation values, the choice
d
neu
2
=
d
dir
should
be favored to have a consistent order in the approximation
of all quantities at the boundaries.
To see that prescribing
d
−
1 points plus some deriva-
tive (resulting in
d
data points) yields the aforementioned
order for the approximation, both a quick argument and
a more formal proof, can be given. As a first qualitative
argument, let’s assume we know the value of a function
f
(
x
) at 0 with an error of order
h
d
, and its first derivative
with error
h
d
−
1
. Then, replacing in its Taylor expansion,
f
(
x
)
≈
f
(0) +
f
′
(0)
h
, it is straightforward to see the er-
ror in both terms is of order
h
d
. More formally, we can
use Lagrange’s mean value theorem to see that there are
d
−
2 intermediate points where the derivatives of the differ-
ence between the function and an interpolating polynomial
vanish. Thus, the derivative of this polynomial is an in-
terpolating polynomial of the derivative of the function
on a set of
d
−
1 points. Then the theorem that gives the
polynomial interpolation error can be applied, and it will
give order
h
d
−
1
. Integrating between 0 and (
d
−
1)
h
gives
the error between the polynomial and the function, giving
an error of order
h
d
. Similar arguments can be written for
the case in which the second derivative is prescribed, or for
the case of Robin boundary conditions.
Finally, it is worth pointing out that in the Dirichlet
case operators
A
dir
and
Q
dir
could be used for any grid,
irrespective of the grid spacing used for generating them.
However, when the second derivative is prescribed, care
should be taken to correctly account for the difference in
the spacings when utilizing
̃
q
neu
2
. Returning now to the
original grid
z
0
, ... ,z
N
−
1
with spacing ∆
z
, and consid-
ering that a
̃
q
neu
2
operator is available and was computed
using a spacing ∆
x
, the value
f
(
z
N
−
1
) can be computed
from
f
′′
(
z
N
−
1
) by simple application of the chain rule, that
is
f
(
z
N
−
1
) =
̃
q
neu
2
f
T
neu
2
,
(43)
with
f
neu
2
=
[
f
(
z
N
−
d
)
, ... ,f
(
z
N
−
2
)
,f
′′
(
z
N
−
1
)(∆
z/
∆
x
)
2
]
.
(44)
7
4.3.
Modified FC-Gram method for Robin boundary condi-
tions
When faced with Robin boundary conditions the func-
tion values
f
(
z
i
) at every interior point are known and the
value of a certain linear combination between the function
and its derivative at one end of the domain is prescribed,
that is
f
′
(
z
N
−
1
) +
λf
(
z
N
−
1
) =
g.
(45)
The first procedure discussed in the preceding subsec-
tion can be also applied for this case. That is, from the
discrete inner product
〈
g
|
h
〉
=
[
g
′
(
x
d
−
1
) +
λg
(
x
d
−
1
)
][
h
′
(
x
d
−
1
) +
λh
(
x
d
−
1
)
]
+
+
d
−
2
∑
i
=0
g
(
x
i
)
h
(
x
i
)
,
(46)
a set of modified Gram polynomials is defined and a related
projector
Q
rob
and blend to zero operator
A
rob
can be
obtained from the following modified Vandermonde matrix
P
rob
=
1
x
0
...
x
d
−
1
0
1
x
1
...
x
d
−
1
1
.
.
.
.
.
.
.
.
.
.
.
.
1
x
d
−
2
...
x
d
−
1
d
−
2
1 1 +
λx
d
−
1
...
(
d
−
1)
x
d
−
2
d
−
1
+
λx
d
−
1
d
−
1
.
(47)
The
C
blend to zero values can then be obtained from
the vector
r
=
[
f
(
x
0
)
,...,f
(
x
d
−
2
)
,f
′
(
x
d
−
1
) +
λf
(
x
d
−
1
)
]
as
A
rob
Q
rob
r
.
It is also possible, and more convenient, to evaluate the
solution value at the endpoint
x
d
−
1
employing the known
values at
x
0
,...,x
d
−
2
and the boundary condition at
x
d
−
1
.
Once again, this can be attained from the relations
P
rob
c
=
r
T
,
(48)
P
dir
c
=
d
T
,
(49)
which has the solution
Q
dir
R
dir
R
−
1
rob
Q
T
rob
r
T
=
d
T
,
(50)
where
R
rob
is the upper triangular matrix resulting from
the
QR
factorization of
P
rob
. The end value
f
(
x
d
−
1
) is
then given by
f
(
x
d
−
1
) =
̃
q
rob
r
T
.
(51)
Here
̃
q
rob
is the last row of
Q
dir
R
dir
R
−
1
rob
Q
T
rob
, that is
( ̃
q
rob
)
i
= (
Q
dir
R
dir
R
−
1
rob
Q
T
rob
)
d
−
1
,i
.
Although more convenient than the method first pro-
posed, this latter algorithm has the limitation that the
matrices
Q
rob
and
R
rob
depend on the specific value of
λ
in the Robin boundary condition, requiring a different high
precision
QR
factorization for each value of
λ
. This would
make this approach prohibitive for the problem considered
in Section 3.2, where a different value of the Robin coupling
parameter
λ
is satisfied by each Fourier mode, and hence
O
(
N
x
N
y
/
4) different operators are required. Changing
the spatial resolution or the domain size would also forbid
reusing operators, as the change in grid spacing also results
in a modification in the value of
λ
. This can be seen from
the chain rule relation
d
f
(
x
)
d
x
+
λf
(
x
) =
g
⇐⇒
d
f
(
x
′
)
d
x
′
+
d
x
d
x
′
λ
︸
︷︷
︸
λ
′
f
(
x
′
) =
d
x
d
x
′
g
︸
︷︷
︸
g
′
.
(52)
Thus, using a different grid spacing requires not only scaling
of the boundary value
g
, as it was in the case of Section 4.2,
but also of the Robin coupling coefficient.
However, the aforementioned limitation can be circum-
vented by proposing the decomposition
P
rob
=
P
neu
+
λ
ˆ
P
dir
,
with
P
neu
the matrix
P
rob
for the special case
λ
= 0 (cor-
responding to a Neumann boundary condition on the first
derivative), and
ˆ
P
dir
the Vandermonde matrix
P
dir
but
replaced with zeros in all but its last row. Substituting
this decomposition in Eqs. (48) and (49) and solving for
c
leads to the equivalent relation
(
Q
neu
R
neu
+
λ
ˆ
P
dir
)
R
−
1
dir
Q
T
dir
f
=
f
rob
,
(53)
where, naturally,
Q
neu
R
neu
is the
QR
factorization of
P
neu
.
Noting that
ˆ
P
dir
can be also represented as
ˆ
Q
dir
R
dir
with
ˆ
Q
the
QR
factorization of
P
dir
but with zeros in all but its
last row, this last expression is reduced to
(
Q
neu
R
neu
R
−
1
dir
Q
T
dir
+
λ
ˆ
I
)
f
=
f
rob
,
(54)
where
ˆ
I
is a matrix whose only non-zero element is
ˆ
I
d
−
1
,d
−
1
=
1. This relation can then be inverted by considering
a
ˆ
I
=
u
T
v
with
u
= (0
,...,
0
,λ
) and
v
= (0
,...,
0
,
1) and
using the Sherman-Morrison formula [40, 41]
(
A
+
uv
T
)
−
1
=
A
−
1
−
A
−
1
uv
T
A
−
1
1 +
v
T
A
−
1
u
,
(55)
which when applied to Eq. (54) and after defining (
̃
q
der
)
i
=
(
Q
dir
R
dir
R
−
1
neu
Q
T
neu
)
d
−
1
,i
, leads to
f
d
−
1
=
(
1
−
a
̃
q
d
−
1
1 +
λ
̃
q
d
−
1
)
̃
q
der
f
T
rob
.
(56)
The obtained value
f
d
−
1
is a
d
-th order approximation to
f
(
x
d
−
1
)
, a fact that can be rigorously demonstrated using
the same arguments provided in Section 4.2.
This allows efficient computation of FC-Gram trans-
forms with Robin boundary conditions with a few pre-
computed coefficients, and without extra costs when do-
main sizes or spatial resolutions are changed. The methods
for a non-periodic function are illustrated in Fig. 2, where
5
th
order (
d
= 5) FC-Gram continuations are applied to
the function
J
0
(25
z
)
e
3
z
−
1 over a grid of
N
= 113 points
spawning the [0
,
1] domain.
8
Algorithm 1:
Schematization of the solver for the perfectly conducting walls scenario.
Starting from the field values at time
t
, (
a
t
,
u
t
,
p
t
,
φ
t
), do:
1:
Obtain intermediate variables
a
∗
t
+∆
t
,
u
∗
t
+∆
t
from the pressureless momentum
equation and the induction equation without the scalar potential.
2:
Apply conditions to the intermediate fields:
a.
homogeneous boundary conditions for
a
∗
t
+∆
t
‖
;
b.
set the mean value of
a
∗
t
+∆
t
⊥
to zero;
c.
non-slip compatible boundary conditions for the velocity.
3:
Solve two Poisson equations to remove the non-solenoidal components in
u
∗
t
+∆
t
,
a
∗
t
+∆
t
:
a.
for the scalar potential
φ
t
+∆
t
with homogeneous Dirichlet boundary conditions;
b.
for the pressure
p
t
+∆
t
with Neumann boundary conditions in order
to cancel the normal velocity at the boundary at projection time.
4:
Perform a solenoidal projection to obtain the fields at the next timestep; that is:
a.
subtract
∇
φ
t
+∆
t
to
a
∗
t
+∆
t
;
b.
subtract
∇
p
t
+∆
t
to
u
∗
t
+∆
t
.
5:
Apply homogeneous boundary conditions to
∂
2
zz
a
t
+∆
t
‖
and
∂
z
a
t
+∆
t
⊥
.
End of iteration
5. A new numerical method for the incompressible
MHD equations with boundaries
We now present a pseudo-spectral method for evolving
the incompressible MHD equations
∂
u
∂t
+ (
u
·
∇
)
u
=
−∇
p
−∇
2
a
×
(
∇
×
a
) +
ν
∇
2
u
,
(57)
∇
2
p
=
−
∇
·
[
(
u
·
∇
)
u
]
,
(58)
∂
a
∂t
=
u
×
(
∇
×
a
) +
η
∇
2
a
−
∇
φ,
(59)
∇
2
φ
=
∇
·
[
u
×
(
∇
×
a
)
]
,
(60)
with the boundary conditions previously discussed. It
should be noted that a Fourier representation is specially
convenient for solving Eqs. (58) and (60), as inverting the
Laplacian in the Fourier domain is attained by simply
inverting a diagonal matrix, contrary to the denser repre-
sentations found in other schemes, such as finite differences
or other spectral methods such as in classic Chebyshev
polynomial decompositions [10].
In all cases, boundary conditions for the velocity field
are periodic in
x
and
y
, and no-normal velocity, no-slip in
z
(i.e.,
u
= 0 in impermeable boundaries).
5.1. Perfectly conducting boundary conditions
We consider the case with two periodic coordinates (
x
and
y
), and a perfect conductor at
z
= 0 and
L
z
. The per-
fect conductor boundary conditions for the vector potential
were discussed in Eqs. (15) to (18), and are repeated here
for convenience:
a
‖
∣
∣
∣
z
=0
,L
z
=
0
,
(61)
∂
2
a
‖
∂z
2
∣
∣
∣
∣
∣
z
=0
,L
z
=
0
,
(62)
∂a
z
∂z
∣
∣
∣
∣
z
=0
,L
z
= 0
.
(63)
As it was explained in Section 3, Eqs. (61) and (62)
are required to represent the physical boundary condition
j
‖
|
z
=0
,L
z
=
0
, which reduces to
∇
2
a
‖
|
z
=0
,L
z
=
0
when the
gauge condition
∇
·
a
= 0 is fixed. Equation (63) is just a
gauge choice whose role is to help in maintaining the error
in the gauge at the boundary small. Moreover, it should be
noted that an error in the condition (
∂
2
zz
a
‖
)
|
z
=0
,L
z
leads
to a non-zero value for
a
‖
|
z
=0
,L
z
, as it can be seen for the
evolution equation for
a
at the boundaries,
∂
a
‖
∂t
∣
∣
∣
∣
z
=0
,L
z
=
η
[
∇
2
a
‖
]
z
=0
,L
z
.
(64)
The latter being a difussion equation, it is hence the case
that any numerical error that might arise in the imposition
of Eq. (62) remains bounded and controlled when the MHD
equations are evolved in time.
To impose Eqs. (61) to (63) we now propose an explicit
time-splitting integration method that imposes boundary
conditions on both the tangential components of
a
‖
and its
second normal derivative, while also utilizing a time-split
scheme for the velocity field as described in [
19
], based on
a method presented in [
42
,
43
]. For simplicity, it will be
described for a first order forward Euler time stepping, but
it can be generalized to any order Runge-Kutta methods in
a way fully analogous to that previously described in [
19
]
(examples in Section 6 were integrated with such a second
order time-splitting Runge-Kutta method).
9