We discuss the application of computer algebra to problems commonly arising
in numerical relativity, such as the derivation of 3+1-splits, manipulation of
evolution equations and automatic code generation. Particular emphasis is put
on working with abstract index tensor quantities as much as possible.Comment: 5 pages, no figures. To appear in the Proceedings of the Spanish
  Relativity Meeting (ERE 2002), Mao, Menorca, Spain, 22-24 Sept 200

Husa, Sascha

Lechner, Christiane

English

arXiv

We discuss the application of computer algebra to problems commonly arising in numerical relativity, such as the derivation of 3+1-splits, manipulation of evolution equations and automatic code generation. Particular emphasis is put on working with abstract index tensor quantities as much as possible

Husa, S.

Lechner, C.

MPG.PuRe

Computer Algebra Applications for Numerical Relativity

arXiv:gr-qc/0301076 v1   21 Jan 2003Computer Algebra Applicationsfor Numerical RelativitySascha Husa1,2 and Christiane Lechner11Max-Planck-Institut fu¨r Gravitationsphysik, Albert-Einstein-Institut, Germany2Departament de F´ısica, Universitat de les Illes Balears, SpainWe discuss the application of computer algebra to problems commonly arisingin numerical relativity, such as the derivation of 3+1-splits, manipulation ofevolution equations and automatic code generation. Particular emphasis is puton working with abstract index tensor quantities as much as possible.Numerical relativists rarely talk about computer algebra (CA) in public, andoutsiders to the field even might get the impression that CA and numerical relativity(NR) indeed represent two antipodal uses of computers to study relativity theory.Quite the opposite is true. In this article we will try to outline some of the possibleapplications of computer algebra for NR, and very briefly present some of our recentwork resulting in the creation of a suite of Mathematica scripts which we have foundextremely useful and which are available on request. While one of our aims is tomake practitioners of the field more aware of the opportunities, another is to helpoutsiders better understand the problems faced in NR. These are often analytical innature and sometimes closer to mainstream mathematical relativity than expected:There is much more to NR than coding up, say, the ADM equations [1]:Lnhab = 2αKab, LnKab = α;ab + α (2KcbKac −KKab − αRab)– or any other particular evolution system one fancies. The perspective of numericalapproximation raises many new questions about the Einstein equations, such as whathappens to the constraints in a free evolution scheme. Mathematical analysis andNR experience have shown that the Einstein equations have to be brought into aform which is suitable for numerical treatment. Considering evolution problems,note that obtaining a well-posed problem is not sufficient. Well posedness doesnot rule out exponential growth which may result from constraint violating modesor a bad gauge. Curing such problems typically requires modifying the equations,and the analysis and coding of different systems of equations. NR thus providesperfect problems for CA, such as: (i) 3+1 or 2+2 decompositions, (ii) modification ofequations by adding constraints, changing variables, etc., (iii) derivation and analysisof associated systems like the constraint propagation system, (iv) linearization aroundexact solutions or (v) the generation of numerical code from systems of equations.12 S. Husa and C. LechnerAll of these tasks can in principle be accomplished by component calculations,as can be carried out quite conveniently and efficiently by computer algebra systemslike GRTensorII [2], which allows to enter expressions in abstract index notation andyields results in component form. This is often what one wants, but thinking aboutderiving and analyzing evolution systems, it is clear that apart from the fact thatthis may result in very large calculations requiring significant time and memory, thismethod is unwieldy and not very intuitive. Rather one would like to keep an abstractindex notation as long as possible, and in particular get results in this notation.The final aim of computer algebra calculations in numerical relativity will usuallybe the generation of code, in C or Fortran say, so let us consider a minimal task list forcode generation. From our pint of view, the problem of generating finite-differenced(or otherwise discretized equations) roughly splits into the following steps:1. Write or find a CA system capable of abstract index tensor calculus.2. Write a package to facilitate 3+1 splits and other calculations in a 3+1 context.3. Derive a system from a 3+1 split or transform a given 3+1 system in some way.At the end of this process the desired system is given in terms of its dependentvariables, their time derivatives and spatial ordinary derivatives.4. Translate the tensor expressions into components.5. Replace the ordinary derivatives by some standard language, e.g. D2[h13].6. Create discretized expressions, e.g.D2[h13] → ( h13(i, j+1, k) - h13(i, j-1, k) ) / dy.7. Wrap this up by code needed to create a full executable programFor the task of coding a simple given system, such as ADM, steps 1 – 3 might beconsidered overkill – it is easy to type in the desired equations in abstract form byhand. For more complicated first order systems or for deriving associated systemslike the constraint propagation system, linearizations or perturbation formalisms,these techniques potentially save a lot of valuable time.As far as code generation is concerned, we decided to generate code as completeCactus [?] thorns. This choice yields an open and documented infrastructure, paral-lelization, clean I/O methods and allows easy interfacing with a growing communitywriting NR Cactus applications. Modifications to interface with other systems withcapabilities similar to those of cactus or with home-brewed code should be ratherstraightforward.For the choice of CA-system we contemplated the use of Maxima, Maple andMathematica. Maxima is a an open source version of Macsyma, and via the itensorpackage supports abstract index calculations. However, to our knowledge itensorComputer Algebra Applications for Numerical Relativity 3is currently not fully functional. Provided development continues, Maxima couldhowever become a very interesting option. Maple and Mathematica are both widelyspread commercial CA systems. Both provide support for component calculations,we find the GRTensorII package for Maple particularly useful. However, as opposedto Mathematica, we are not aware of any functional abstract index tensor pack-age for Maple. We speculate that this is rooted in Mathematica’s superior intrinsicsupport for pattern matching. This seems quite essential for tensor manipulations,e.g. Tab and Tcd are not the same expression but still are equivalent mathemati-cal objects, which can easily be identified with pattern matching techniques. Wehave worked with two different abstract index packages for Mathematica: the freelyavailable Ricci and the commercial MathTensor. Despite our feeling that the overalldesign of MathTensor is less clean, and despite the fact that we have found severalbugs and inconsistencies in MathTensor, we still eventually selected it as the basisfor our work. The main reasons were a somewhat more extensive functionality, itssupport for both abstract and component calculations, better checking for errors(such as inconsistencies with indices), and it’s simpler representation of tensors asplain functions of indices: h[la, lb]→ hab, CD[Metricg[la, lb], lc]→ 0. Thisstraightforward syntax is less error-prone than Ricci’s corresponding h[L[a], L[b]]or h[L[a], L[b]] [L[c]] for a covariant derivative. Given MathTensor’s immensevalue, problems and significant cost, it would be very attractive to have available anopen source alternative with similar functionality.The basic set of functions needed for 3+1 decompositions, as well as some generaltensor manipulation functionality, has been implemented in a Mathematica notebookDecompose 3+1 Tools.nb. Our strategy to do 3+1 decompositions using MathTen-sor, can be outlined as follows• Define tensors to be labeled spatial; define a vector (n, the unit normal) to belabeled timelike using functions to generate and manage lists of hypersurface-tensors, e.g. for the ADM equations the calls would look something like:DefineSpatialVector[Shift]; DefineTimelikeVector[n, t]DefineSpatial2Tensors[h, K, 1];• Instruct MathTensor about projection rules such as naT...a → 0, naLnT...a... →0 for T spacelike. Such rules are defined for all spatial tensors by callingDefineHypersurfaceOrthogonalityRules[h, n], where h is the metric in-duced on the hypersurface and n its unit normal.• The function DefineFundamentalFormsRules[h, K, a, n, Dh] attaches namesto the geometrical objects (n is the unit normal, Dh the induced covariantderivative, ab = nc∇cnb) and defines the decomposition of the 4 metric gab =hab + nanb and the definition of the extrinsic curvature ∇anb = δ Kab + naabwith  and δ global variables to serve different sign conventions. Many associ-ated rules get defined automatically.4 S. Husa and C. Lechner• Define additional rules that are particular to the problem treated. Examples forsuch additional rules would be the splitting of the Maxwell tensor into electricand magnetic fields in Maxwell theory. Rules for standard formulas like theGauss-Codazzi relations or the split of the unit normal into lapse and shift aredefined by calling high-level functions.• Compute all independent projections of the 4-dimensional field equations (withrespect to h, n) and use the projection rules defined above, afterwards switchto manifest three dimensional form (e.g. set hji = δji ).• Compute equations for first order variables (e.g. for the Christoffel symbols).We have defined functions to compute Lie- and covariant derivatives for tensor densi-ties in terms of ordinary derivatives, which is not directly supported by MathTensor,e.g. this code defines a rule to deal with Lie derivatives of tensor densities by addingthe appropriate correction term:RuleUnique[densityLieDRule, LieD[T_[xa__], v_], (LieD[T[xa], v])+ densityWeight[T] T[xa] OD[v[uc], lc], MemberQ[densityList, T]]RuleUnique is a MathTensor command to create rules which respect dummy indicesand densityList is a list of all tensor density objects.Mathematica notebooks containing examples for 3+1 decompositions, startingfrom 4-dimensional equations up the the generation of code for evolution systemand evaluation of constraints have been worked out for the Maxwell equations, theADM equations [1], the conformal field equations [1] and the BSSN equations [1].Notebooks working out the constraint propagation system have been developed forMaxwell and the conformal field equations. Treatment of other systems should bestraightforward following these examples.The key to generate code is to generate the lists of independent tensor compo-nents and component equations. Sums over indices are expanded with MathTensor’sMakeSum command. Assigning names to these variables as they should appear inthe code (e.g. {h11, h12, h13, . . .}) is straightforward with Mathematica’s patternmatching techniques. Ordinary derivatives to a standard syntax which can easily beexpanded to finite difference expressions with Macros. Here we use rules such asDiffCompsRule = OD[T_, n_?IntegerQ] :>ToExpression["D"<>ToString[Abs[n]]] @ T;which would yield OD[h32, -2]→ D2h32. Simplifications of the resulting componentexpressions are obtained with Mathematica’s Collect function. Along these lines wedeveloped a function to generate Cactus evolution thorns and similar functions togenerate thorns that evaluate constraints or any user-define geometric quantities orset gauge and initial data from a 4-metric (e.g. an exact solution). Special care hasComputer Algebra Applications for Numerical Relativity 5been taken to generate nicely formatted human readable code and to not assume aparticular system of equations or set of variables.Summing up, we have tried to promote and discuss the use of computer algebrafor NR. CA makes it easier to focus on algorithms, detached from a particular sys-tem of equations. Stressing a more abstract point of view is not only mathematicallymore appealing but also increases flexibility, which benefits scientific productivity.We have very briefly outlined the content of a set of Mathematica scripts which wedeveloped for our use. These scripts consist of functions to manage the bookkeepingof hypersurface-related quantities, the definition of associated rules within MathT-ensor, addition of higher-level tensor and component manipulation functions, codegeneration scripts, and the development of a set of examples to be used as templatesto deal with formulations we have not covered. We have not produced software inthe sense of user-friendly, well documented programs with online help and extensiveerror-checking. Things are not as well automated as they could be – e.g. to dealwith a new system one would essentially follow our template notebooks instead ofcalling just a few high-level functions. Some of our current code is hardwired to 3+1decompositions, as opposed to 2+2 etc. – it would be interesting to generalize ourcode and techniques in this respect, or even to higher dimensions. Despite its defi-ciencies, we consider our work potentially useful for others, and it is freely availableon request. What goes much beyond our scope and resources is to build a reasonablywell documented free community tool for computer algebra in the context of NR,e.g. comparable to the Cactus computational toolkit – it could however lead to acomparable increase in productivity! Needless to say, we have found it invaluable toaccompany our computer calculations by “unplugged” manual work.Acknowledgments: We thank the members of AEI’s Astrophysical NumericalRelativity group for stimulating discussions on the topic, S.H. has benefited fromreading Mathematica code written by Bernd Bru¨gmann and Peter Hu¨bner.References[1] For an overview of different formulations used in NR and for references see e.g.L. Lehner, Class. Quant. Grav. 18, R25 (2001)[2] P. Musgrave, D. Pollney and K. Lake, Fields Institute Comm. 15 313 (1996)(http://grtensor.phy.queensu.ca/)[3] G. Allen et. al., Cluster Computing, 4 179 (2001); B. Talbot, S. Zhou, G. Higgins,NASA report, http://ct.gsfc.nasa.gov/esmf tasc/Files/Cactus b.html[4] J. Lee, Ricci, a mathematica package for doing tensor calculations in differentialgeometry, http://www.math.washington.edu/∼lee[5] L. Parker and S. Christensen, MathTensor: A System for Doing Tensor Analysisby Computer, Addison-Wesley 1994

http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.339.6573

Computer Algebra Applications for Numerical Relativity

Abstract

Similar works

Full text

Available Versions

MPG.PuRe

MPG.PuRe