Geoscience Reference
In-Depth Information
continuation we describe are essentially the same as those
that would be applied to low-dimensional dynamical sys-
tems, except that any dense linear problems that arise must
be solved by iterative methods because direct methods are
prohibitively memory and time intensive.
For the application of continuation methods to a
branch of steady solutions, we could proceed, as we did
in Section 2.3, by setting the time derivative equal to zero.
Instead, we elect to compute steady solutions as fixed
points of the flow over some finite time
t
. This makes
the extension to periodic orbits straightforward and, as
explained below, yields linear problems directly amenable
to iterative solvers.
We define the map
to be the flow of the discretized
system, i.e., the solution
X
(t)
of the discretized model
equations with initial condition
X
(
0
)
, for fixed values of
the continuation parameter
α
:
tangent vector
t
i
. Thus, instead of solving (2.22) alone, we
solve the extended system:
(
X
i
+1
,
t
,
α
i
+1
)
−
X
i
+1
=0,
(2.24)
X
i
+1
)
α
i
+1
)t
iα
=0.
(
X
i
+1
−
·
t
i
X
+
(α
i
+1
−
The second subscript attributed to the tangent vector
denotes its component. The linear system that is to be
solvedatthe
k
thNewton-Raphsoniterationtakestheform
∇
X
X
−
I
∇
α
t
i
X
t
iα
α
(
X
k
−
1
X
k
−
1
i
+1
i
+1
,
t
,
α
i
+1
)
−
=
−
(2.25)
(
X
k
−
1
i
+1
X
i
+1
,
α
k
−
1
α
i
+1
)
t
t
i
−
−
i
+1
and leads to the
k
th approximation to the
(i
+1
)
(i+1)st solution
on the branch:
(
X
i
+1
,
α
i
+1
)
=
(
X
k
−
1
i
+1
,
α
k
−
1
d
dt
(
X
(
0
)
,
t
,
α)
=
F
((
X
(
0
)
,
t
,
α))
,
i
+1
)
+
(
X
,
α)
.
This linear system is large and dense and can, in gen-
eral, not be solved by direct methods. Instead, we can use
an iterative solver like GMRES [
Saad and Schultz
, 1986].
This method constructs an approximate solution to the
linear systems based on a sequence of matrix-vector prod-
ucts. Each product can be computed from an integration
of the linearized equations. Thus, each GMRES iteration
requires one integration of the discretized model, along
with its linearization, over time
t
. The number of iterations
necessary to obtain a solution with sufficient accuracy
to ensure convergence of the Newton-Raphson iteration
depends strongly on the spectrum of the matrix. Gener-
ally speaking, the more clustered the spectrum, the faster
the convergence of GMRES. Because of the dissipative
nature of the discretized model, the spectrum will be more
clustered for longer integration time
t
, since the majority
of the eigenvalues of
(2.21)
where
X
(t)
is a vector formed from the discretized veloc-
ity vector
u
and the temperature
T
,
F
corresponds to the
discretization of equations (2.4)-(2.8), and the continua-
tion parameter
α
could be, for instance, the rotation rate
or the differential heating
T
. The stationarity condition
canthenbewrittenas:
−
(
X
,
t
,
α)
X
=
0
(2.22)
for arbitrary time
t
.
Generically, the solutions lie on curves, i.e., the solu-
tion branches. To proceed, it is necessary to know an
initial point
(
X
0
,
α
0
)
on the solution branch. This can be
a stable fixed point computed by numerical integration
or the discretization of a theoretically known solution.
Given a known point
(
X
i
,
α
i
)
along a branch, the contin-
uation method is composed of two parts: (1) a prediction
step, in which an initial guess
(
X
i
+1
,
α
i
+1
)
of a new point
(
X
i
+1
,
α
i
+1
)
along the branch is formed, and (2) a correc-
tion step, in which the initial guess is iteratively refined
until it is sufficiently close to an exact solution.
The prediction
(
X
i
+1
,
α
i
+1
)
is generated from the previ-
ous (known) point
(
X
i
,
α
i
)
and from an approximation of
the tangent vector to the branch,
t
i
. Specifically, we choose
(
X
i
+1
,
α
i
+1
)
=
(
X
i
,
α
i
)
+
t
i
,
correspond to strongly damped
perturbations and accumulate around zero. Consequently,
the choice of the integration time is a trade-off between
the time taken by each GMRES iteration on one hand and
the required number of iterations on the other.
The algorithm for tracking periodic solutions is largely
the same as that for stationary states. The main difference
is that the integration time is now the period
τ
of the orbit,
which is unique but unknown and therefore becomes a
new variable of the problem. An additional scalar equa-
tion is needed for closure. Given a solution
(
X
i
,
α
i
,
τ
i
)
and
a prediction
(
X
i
+1
,
α
i
+1
,
τ
i
+1
)
for the next solution, the
equations to solve are
∇
(2.23)
where
is the step size, which approximates the increment
of the arc length along the computed solution branch.
The correction step is based on Newton-Raphson
iteration. However, the linear equation to be solved in each
Newton-Raphson iteration is underdetermined, reflecting
the fact that the next fixed point is not uniquely defined
but lies on a curve. To render the next fixed point unique,
we seek a solution in the hyperplane orthogonal to the
(
X
i
+1
,
τ
i
+1
,
α
i
+1
)
−
X
i
+1
=0,
X
i
+1
)
α
i
+1
)t
iα
(
X
i
+1
−
·
t
i
X
+
(α
i
+1
−
τ
i
+1
)t
iτ
= 0,
+
(τ
i
+1
−
(2.26)
f (
X
i
+1
)
=0,
