Tensor calculus with open-source software: the SageManifolds project

The SageManifolds project aims at extending the mathematics software system Sage towards differential geometry and tensor calculus. Like Sage, SageManifolds is free, open-source and is based on the Python programming language. We discuss here some details of the implementation, which relies on Sage's parent/element framework, and present a concrete example of use.


Introduction
Computer algebra for general relativity (GR) has a long history, which started almost as soon as computer algebra itself in the 1960s. The first GR program was GEOM, written by J.G. Fletcher in 1965 [1]. Its main capability was to compute the Riemann tensor of a given metric. In 1969, R.A. d'Inverno developed ALAM (for Atlas Lisp Algebraic Manipulator ) and used it to compute the Riemann and Ricci tensors of the Bondi metric. According to [2], the original calculations took Bondi and collaborators 6 months to finish, while the computation with ALAM took 4 minutes and yielded the discovery of 6 errors in the original paper. Since then, numerous packages have been developed: the reader is referred to [3] for a review of computer algebra systems for GR prior to 2002, and to [4] for a more recent review focused on tensor calculus. It is also worth to point out the extensive list of tensor calculus packages maintained by J. M. Martin-Garcia at [5].
All applications listed in the second category are free software. In the first category, xAct and Ricci are also free software, but they require a proprietary product, the source code of which is closed (Mathematica).

The SageManifolds project 4.1. Aim and scope
Sage is well developed in many areas of mathematics but very little exists for differential geometry and tensor calculus. One may mention differential forms defined on a fixed coordinate patch, implemented by J. Vankerschaver [20], and the 2-dimensional parametrized surfaces of the 3dimensional Euclidean space recently added by M. Malakhaltsev and J. Vankerschaver [21]. The aim of SageManifolds [22] is to introduce smooth manifolds and tensor fields in Sage, with the following requirements: (i) one should be able to introduce various coordinate charts on a manifold, with the relevant transition maps; (ii) tensor fields must be manipulated as such and not through their components with respect to a specific (possibly coordinate) vector frame.
Concretely, the project amounts to creating new Python classes, such as Manifold, Chart, TensorField or Metric, to implement them within the parent/element pattern and to code mathematical operations as class methods. For instance the class Manifold, devoted to real smooth manifolds, is a parent class, i.e. it inherits from Sage's class Parent. On the other hand, the class devoted to manifold points, ManifoldPoint, is an element class and therefore inherits from Sage's class Element. This is illustrated by the inheritance diagram of Fig. 1  diagram, each class at the base of some arrow is a subclass (also called derived class) of the class at the arrowhead. Note however that the actual type of a parent is a dynamically generated class taking into account the mathematical category to which it belongs. For instance, the actual type of a smooth manifold is not Manifold, but a subclass of it named Manifold with category, reflecting the fact that Manifold is declared in the category of Sets 2 . Note also that the class Manifold inherits from Sage's class UniqueRepresentation, which ensures that there is a unique manifold instance for a given dimension and given name.

Implementation of charts
Given a smooth manifold M of dimension n, a coordinate chart on some open subset U ⊂ M is implemented in SageManifolds via the class Chart, whose main data is a n-tuple of Sage symbolic variables (x 1 , . . . , x n ), each of them representing a coordinate. In general, more than one (regular) chart is required to cover the entire manifold. For instance, at least 2 charts are necessary for the n-dimensional sphere S n (n ≥ 1) and the torus T 2 and 3 charts for the real projective plane RP 2 (see Fig. 6 below). Accordingly, SageManifolds allows for an arbitrary number of charts. To fully specify the manifold, one shall also provide the transition maps (changes of coordinates) on overlapping chart domains (SageManifolds class CoordChange).

Implementation of scalar fields
A scalar field on manifold M is a smooth mapping  where U is some open subset of M. A scalar field has different coordinate representations F , F , etc. in different charts X,X, etc. defined on U : These representations are stored in some attribute of the class ScalarField, namely a Python dictionary 3 whose keys are the various charts defined on U : f. express = X : F,X :F , . . . .
Each representation F is an instance of the class FunctionChart, which resembles Sage native symbolic functions, but involves automatic simplifications in all arithmetic operations. Given an open subset U ⊂ M, the set C ∞ (U ) of scalar fields defined on U has naturally the structure of a commutative algebra over R: it is clearly a vector space over R and it is endowed with a commutative ring structure by pointwise multiplication: The algebra C ∞ (U ) is implemented in SageManifolds via the parent class ScalarFieldAlgebra, in the category CommutativeAlgebras.
The corresponding element class is of course ScalarField (cf. Fig. 2).

Modules and free modules
Given an open subset U ⊂ M, the set X (U ) of all smooth vector fields defined on U has naturally the structure of a module over the algebra C ∞ (U ). Let us recall that a module is similar to a vector space, except that it is based on a ring (here C ∞ (U )) instead of a field (usually R or C in physical applications). Of course, every vector space is a module, since every field is a ring. There is an important difference though: every vector space has a basis (as a consequence of the axiom of choice), while a module does not necessarily have any. When it possesses one, it is called a free module. Moreover, if the module's base ring is commutative, it can be shown that all bases have the same cardinality, which is called the rank of the module (for vector spaces, which are free modules, the word dimension is used instead of rank ).
If X (U ) is a free module (examples are provided in Sec. 4.5 below), a basis of it is nothing but a vector frame (e a ) 1≤a≤n on U (often called a tetrad in the context of 4-dimensional GR): The rank of X (U ) is thus n, i.e. the manifold's dimension 4 . At any point p ∈ U , Eq. (5) gives birth to an identity in the tangent vector space T p M: which means that the set (e a (p)) 1≤a≤n is a basis of T p M. Note that if U is covered by a chart (x a ) 1≤a≤n , then (∂/∂x a ) 1≤a≤n is a vector frame on U , usually called coordinate frame or natural basis. Note also that, being a vector space over R, the tangent space T p M represents another kind of free module which occurs naturally in the current context. It turns out that so far only free modules with a distinguished basis were implemented in Sage. This means that, given a free module M of rank n, all calculations refer to a single basis of M . This amounts to identifying M with R n , where R is the ring over which M is based. This is unfortunately not sufficient for dealing with smooth manifolds in a coordinate-independent way. For instance, there is no canonical isomorphism between T p M and R n when no coordinate system is privileged in the neighborhood of p. Therefore we have started a pure algebraic part of SageManifolds to implement generic free modules, with an arbitrary number of bases, none of them being distinguished. This resulted in (i) the parent class FiniteRankFreeModule, within Sage's category Modules, and (ii) the element class FiniteRankFreeModuleElement. Then both classes VectorFieldFreeModule (for X (U ), when it is a free module) and TangentSpace (for T p M) inherit from FiniteRankFreeModule (see Fig. 3).

Implementation of vector fields
Ultimately, in SageManifolds, vector fields are described by their components with respect to various vector frames, according to Eq. (5), but without any vector frame being privileged, leaving the freedom to select one to the user, as well as to change coordinates. A key point is that not every manifold admits a global vector frame. A manifold M, or more generally an open subset U ⊂ M, that admits a global vector frame is called parallelizable. Equivalently, M is parallelizable if, and only if, X (M) is a free module. In terms of tangent bundles, parallelizable manifolds are those for which the tangent bundle is trivial: Examples of parallelizable manifolds are [23] • the Cartesian space R n for n = 1, 2, . . ., • the sphere S 3 SU(2), as any Lie group, • the sphere S 7 , • any orientable 3-manifold (Steenrod theorem [24]).  • the sphere S 2 (as a consequence of the hairy ball theorem), as well as any sphere S n with n ∈ {1, 3, 7}, • the real projective plane RP 2 .
Actually, "most" manifolds are non-parallelizable. As noticed above, if a manifold is covered by a single chart, it is parallelizable (the prototype being R n ). But the reverse is not true: S 1 and T 2 are parallelizable and require at least two charts to cover them.
If the manifold M is not parallelizable, we assume that it can be covered by a finite number In particular, this holds if M is compact, for any compact manifold admits a finite atlas. We then consider the restrictions of vector fields to the U i 's. For each i, X (U i ) is a free module of rank n = dim M and is implemented in SageManifolds as an instance of VectorFieldFreeModule (cf. Sec 4.6. Implementation of tensor fields The implementation of tensor fields in SageManifolds follows the strategy adopted for vector fields. Consider for instance a tensor field T of type (1,1) on the manifold M. It can which defines T a b , is meaningful only when a vector frame (e a ) exists 5 . Therefore, one first decomposes the tensor field T into its restrictions T | U i on parallelizable open subsets of M, U i (1 ≤ i ≤ N ) and then considers the components on various vector frames on each subset U i . For each vector frame (e a ), the set of components (T a b ) is stored in a devoted class (named Components), which takes into account all the tensor monoterm symmetries: only non-redundant components are stored, the other ones being deduced by (anti)symmetry. This is illustrated in Fig. 5, which depicts the internal storage of tensor fields in SageManifolds. Note that each component T a b is a scalar field on U i , according to the formula Accordingly, the penultimate level of Fig. 5 corresponds to the scalar field storage, as described by (3). The last level is constituted by Sage's symbolic expressions (class Expression).

Parallelization
To improve the reactivity of SageManifolds and take advantage of multicore processors, some tensorial operations are performed by parallel processes. The parallelization is implemented by means of the Python library multiprocessing, via the built-in Sage decorator @parallel. Using it permits to define a function that is run on different sub-processes. If n processes are used, given a function and a list of arguments for it, any process will call the function with an element of the list, one at time, spanning all the list.
Currently 6 , the parallelized operations are tensor algebra, tensor contractions, computation of the connection coefficient and computation of Riemann tensor.
The parallelization of an operation is achieved by first creating a function which computes the required operation on a subset of the components of a tensor; second, by creating a list of 2n (twice the number of used processes) arguments for this function. Then applying this function to the input list, the calculation is performed in parallel. At the end of the computation a fourth phase is needed to retrieve the results. The choice to divide the work in 2n is a compromise between the load balancing and the cost of creating multiple processes. The number of processors to be used in the parallelization can be controlled by the user.

SageManifolds at work: Kerr spacetime and Simon-Mars tensor
We give hereafter a short illustration of SageManifolds focused on tensor calculus in 4-dimensional GR. Another example, to be found at [25], is based on the manifold S 2 and focuses more on the use of multiple charts and on the treatment of non-parallelizable manifolds. Yet another example illustrates some graphical capabilities of SageManifolds: Figure 6 shows the famous immersion of the real projective plane RP 2 into the Euclidean space R 3 known as the Boy surface. This figure has been obtained by means of the method plot() applied to three coordinate charts covering RP 2 , the definition of which is related to the interpretation of RP 2 as the set of straight lines ∆ through the origin of R 3 : (i) in red, the chart X 1 covering the open subset of RP 2 defined by all lines ∆ that are not parallel to the plane z = 0, the coordinates of X 1 being the coordinates (x, y) of the intersection of the considered line ∆ with the plane z = 1; (ii) in green, the chart X 2 covering the open subset defined by all lines ∆ that are not parallel to the plane x = 0, the coordinates of X 2 being the coordinates (y, z) of intersection with the plane x = 1; (iii) in blue, the chart X 3 covering the open subset defined by all lines ∆ that are not parallel to the plane y = 0, the coordinates of X 2 being the coordinates (z, x) of intersection with the plane y = 1. Figure 6 actually shows the coordinate grids of these three charts through the Apéry map [26], which realizes an immersion of RP 2 into R 3 . This example, as many others, can be found at [25].

function(arguments)
In a functional language, this would correspond to result = function(object,arguments). For instance, the Riemann tensor of a metric g is obtained as riem = g.riemann() (in this case, there is no extra argument, hence the empty parentheses). With this in mind, let us proceed with the computation by means of SageManifolds. In the text below, the blue color denotes the outputs as they appear in the Sage notebook (note that all outputs are automatically L A T E Xformatted by Sage).
The first step is to declare the Kerr spacetime (or more precisely the part of the Kerr spacetime covered by Boyer-Lindquist coordinates) as a 4-dimensional manifold:
Let us conclude by stating that, in the very spirit of free software, anybody interested in contributing to the project is very welcome!