Philip W. Livermore

Spectral radial basis functions for full sphere computations

The singularity of cylindrical or spherical coordinate systems at the origin imposes certain regularity conditions on the spectral expansion of any infinitely differentiable function. There are two efficient choices of a set of radial basis functions suitable for discretising the solution of a partial differential equation posed in either such geometry. One choice is methods based on standard Chebyshev polynomials; although these may be efficiently computed using fast transforms, differentiability to all orders of the obtained solution at the origin is not guaranteed. The second is the so-called one-sided Jacobi polynomials that explicitly satisfy the required behavioural conditions. In this paper, we compare these two approaches in their accuracy, differentiability and computational speed. We find that the most accurate and concise representation is in terms of one-sided Jacobi polynomials. However, due to the lack of a competitive fast transform, Chebyshev methods may be a better choice for some computationally intensive timestepping problems and indeed will yield sufficiently (although not infinitely) differentiable solutions provided they are adequately converged.

Galerkin orthogonal polynomials

Abstract The Galerkin method offers a powerful tool in the solution of differential equations and function approximation on the real interval [−1, 1]. By expanding the unknown function in appropriately chosen global basis functions, each of which explicitly satisfies the given boundary conditions, in general this scheme converges exponentially fast and almost always supplies the most terse representation of a smooth solution. To date, typical schemes have been defined in terms of a linear combination of two Jacobi polynomials. However, the resulting functions do not inherit the expedient properties of the Jacobi polynomials themselves and the basis set will not only be non-orthogonal but may, in fact, be poorly conditioned. Using a Gram-Schmidt procedure, it is possible to construct, in an incremental fashion, polynomial basis sets that not only satisfy any linear homogeneous boundary conditions but are also orthogonal with respect to the general weighting function ( 1 - x ) α ( 1 + x ) β . However, as it stands, this method is not only cumbersome but does not provide the structure for general index n of the functions and obscures their dependence on the parameters ( α , β ) . In this paper, it is shown that each of these Galerkin basis functions, as calculated by the Gram-Schmidt procedure, may be written as a linear combination of a small number of Jacobi polynomials with coefficients that can be determined. Moreover, this terse analytic representation reveals that, for large index, the basis functions behave asymptotically like the single Jacobi polynomial P n ( α , β ) ( x ) . This new result shows that such Galerkin bases not only retain exponential convergence but expedient function-fitting properties too, in much the same way as the Jacobi polynomials themselves. This powerful methodology of constructing Galerkin basis sets is illustrated by many examples, and it is shown how the results extend to polar geometries. In exploring more generalised definitions of orthogonality involving derivatives, we discuss how a large class of differential operators may be discretised by Galerkin schemes and represented in a sparse fashion by the inverse of band-limited matrices.

Successive elimination of shear layers by a hierarchy of constraints in inviscid spherical-shell flows

In a rotating spherical shell, an inviscid inertia-free flow driven by an arbitrary body force will have cylindrical components that are either discontinuous across, or singular on, the tangent cylinder, the cylinder tangent to the inner core and parallel to the rotation axis. We investigate this problem analytically, and show that there is an infinite hierarchy of constraints on this body force which, if satisfied, sequentially remove discontinuities or singularities in flow derivatives of progressively higher order. By splitting the solution into its equatorial symmetry classes, we are able to provide analytic expressions for the constraints and demonstrate certain inter-relations between them. We show numerically that viscosity smoothes any singularity in the azimuthal flow component into a shear layer, comprising inner and outer layers, either side of the tangent cylinder, of width O(E 2/7) and O(E 1/4), respectively, where E is the Ekman number. The shear appears to scale as O(E −1/3) in the equatorially symmetric case, although in a more complex fashion when considering equatorial antisymmetry, and attains a maximum value in either the inner or outer sublayers depending on equatorial symmetry. In the low-viscosity magnetohydrodynamic system of the Earths core, magnetic tension within the fluid resists discontinuities in the flow and may dynamically adjust the body force in order that a moderate number of the constraints are satisfied. We speculate that it is violations of these constraints that excites torsional oscillations, magnetohydrodynamic waves that are observed to emanate from the tangent cylinder.

On magnetic energy instability in spherical stationary flows

Eigenmode or linear analysis of the induction equation yields magnetic field structures which grow or decay exponentially in time, under the influence of a prescribed flow geometry. However, due to the non-self-adjointness of the problem, interference between the non-orthogonal modes can lead to subcritical transient growth, the onset of which is described by critical magnetic energy stability. Using a variational principle, the energy method is detailed and for several different spherical flows, calculations are presented which indicate robust results, in contrast to linear analysis which depends critically on the precise definition of the flow used. All flows studied indicate an apparent linear asymptotic dependence of energetic instability on Rm, the magnetic Reynolds number, which is reached when Rm is O(1000). For the flows studied, we find improved lower bounds on Rm for energetic instability of between 5 and 10 times, compared with the analytic bound of Proctor.

The structure of Taylor's constraint in three dimensions

In a 1963 edition of Proc. R. Soc. A, J. B. Taylor (Taylor 1963 Proc. R. Soc. A 9, 274–283) proved a necessary condition for dynamo action in a rapidly rotating electrically conducting fluid in which viscosity and inertia are negligible. He demonstrated that the azimuthal component of the Lorentz force must have zero average over any geostrophic contour (i.e. a fluid cylinder coaxial with the rotation axis). The resulting dynamical balance, termed a Taylor state, is believed to hold in the Earths core, hence placing constraints on the class of permissible fields in the geodynamo. Such states have proven difficult to realize, apart from highly restricted examples. In particular, it has not yet been shown how to enforce the Taylor condition exactly in a general way, seeming to require an infinite number of constraints. In this work, we derive the analytic form for the averaged azimuthal component of the Lorentz force in three dimensions after expanding the magnetic field in a truncated spherical harmonic basis chosen to be regular at the origin. As the result is proportional to a polynomial of modest degree (simply related to the order of the spectral expansion), it can be made to vanish identically on every geostrophic contour by simply equating each of its coefficients to zero. We extend the discussion to allow for the presence of an inner core, which partitions the geostrophic contours into three distinct regions.

Electromagnetically driven westward drift and inner-core superrotation in Earth's core.

Significance Seismic probing of the earth’s deep interior has shown that the inner core, the solid core of our planet, rotates slightly faster (i.e., eastward) than the rest of the earth. Quite independently, observations of the geomagnetic field provide evidence of westward-drifting features at the edge of the liquid outer core. This paper describes a computer model that suggests that the geomagnetic field itself may provide a link between them: The associated electromagnetic torque currently is westward in the outermost outer core, whereas an equal and opposite torque is applied to the inner core. Decadal changes in the geomagnetic field may cause fluctuations in both these effects, consistent with recent observations of a quasi-oscillatory inner-core rotation rate. A 3D numerical model of the earth’s core with a viscosity two orders of magnitude lower than the state of the art suggests a link between the observed westward drift of the magnetic field and superrotation of the inner core. In our model, the axial electromagnetic torque has a dominant influence only at the surface and in the deepest reaches of the core, where it respectively drives a broad westward flow rising to an axisymmetric equatorial jet and imparts an eastward-directed torque on the solid inner core. Subtle changes in the structure of the internal magnetic field may alter not just the magnitude but the direction of these torques. This not only suggests that the quasi-oscillatory nature of inner-core superrotation [Tkalčić H, Young M, Bodin T, Ngo S, Sambridge M (2013) The shuffling rotation of the earth’s inner core revealed by earthquake doublets. Nat Geosci 6:497–502.] may be driven by decadal changes in the magnetic field, but further that historical periods in which the field exhibited eastward drift were contemporaneous with a westward inner-core rotation. The model further indicates a strong internal shear layer on the tangent cylinder that may be a source of torsional waves inside the core.

An implementation of the exponential time differencing scheme to the magnetohydrodynamic equations in a spherical shell

Over the last decade there has been renewed interest in applying exponential time differencing (ETD) time stepping schemes to the solution of stiff systems. In this paper, we present an implementation of such a scheme to the fully spectral solution of the incompressible magnetohydrodynamic equations in a spherical shell. One problem associated with ETD schemes is the accurate calculation of the necessary matrices; we implement and discuss in detail a variety of different methods including direct computation, contour integration, spectral expansions and recurrence relations. We compare the accuracy of six different second-order methods in determining the evolution of a three-dimensional magnetic field under the action of a prescribed time-dependent flow of electrically conducting fluid, and find that for the timestep restriction imposed by the nonlinear terms, ETD methods are no more accurate than linearly implicit methods which have the significant advantage of being easier to implement. However, ETD methods are more readily extendable than those which are linearly implicit and will become much more advantageous at higher order.

Quasi-L p norm orthogonal Galerkin expansions in sums of Jacobi polynomials: Orthogonal expansions

In the study of differential equations on [ − 1,1] subject to linear homogeneous boundary conditions of finite order, it is often expedient to represent the solution in a Galerkin expansion, that is, as a sum of basis functions, each of which satisfies the given boundary conditions. In order that the functions be maximally distinct, one can use the Gram-Schmidt method to generate a set orthogonal with respect to a particular weight function. Here we consider all such sets associated with the Jacobi weight function, w(x) = (1 − x)α(1 + x)β. However, this procedure is not only cumbersome for sets of large degree, but does not provide any intrinsic means to characterize the functions that result. We show here that each basis function can be written as the sum of a small number of Jacobi polynomials, whose coefficients are found by imposing the boundary conditions and orthogonality to the first few basis functions only. That orthogonality of the entire set follows—a property we term “auto-orthogonality”—is remarkable. Additionally, these basis functions are shown to behave asymptotically like individual Jacobi polynomials and share many of the latter’s useful properties. Of particular note is that these basis sets retain the exponential convergence characteristic of Jacobi expansions for expansion of an arbitrary function satisfying the boundary conditions imposed. Further, the associated error is asymptotically minimized in an Lp(α) norm given the appropriate choice of α = β. The rich algebraic structure underlying these properties remains partially obscured by the rather difficult form of the non-standard weighted integrals of Jacobi polynomials upon which our analysis rests. Nevertheless, we are able to prove most of these results in specific cases and certain of the results in the general case. However a proof that such expansions can satisfy linear boundary conditions of arbitrary order and form appears extremely difficult.

Transient magnetic energy growth in spherical stationary flows

The non-normality associated with the induction equation may lead to subcritical growth of magnetic field structures even if all linear eigenmodes decay. We compute the magnitude of these transient effects for a collection of predominantly axisymmetric stationary spherical flows and find without exception the dominance of axisymmetric field growth above all other symmetries. The transient growth is robust under small flow perturbations and can be understood by simple physical mechanisms: either field line shearing or stretching. Magnetic energy amplification of is possible at magnetic Reynolds numbers of , and such effects could therefore lead the system from a principally non-magnetic state into one where the magnetic field plays a significant role in the dynamics.

The role of helicity and stretching in forced kinematic dynamos in a spherical shell

Considerations of mean-field theory suggest that small-scale helical flows are an effective means of generating large-scale (mean) magnetic fields, whereas fast dynamo considerations reveal the importance of Lagrangian chaos in the flow for generating small-scale magnetic fields in the limit of high magnetic Reynolds number. We explore these ideas further by considering the kinematic magnetic fields generated by three forced steady flows in a spherical shell that differ both in their helicity and in their stretching properties. The full magnetic induction equation is solved numerically, with no a priori assumptions about the nature of the generated magnetic field. There are two surprising aspects to our results. One is that the most significant mean field is generated by a flow with zero net helicity; the other is that the flow with the “best” stretching properties turns out to be the most inefficient dynamo. Our results, therefore, suggest that it may not be possible to determine the nature of a kinemati...