Finite Element Analysis & General Relativity

Many years ago I studied the applicability of an unusual computer architecture to computational problems in hydrodynamics. I stumbled onto some insights that might bear on computations in general relativity that might be useful or at least amusing.

They year was 1961. The machine was the Harvest (IBM 7950). Harvest was adapted to processing arrays of data. It was efficient at addition, subtraction, table lookup, and a variety of boolean operations. It was very flexible regarding the number of bits in a datum and provided nearly constant bit rate performance thus rewarding short operands.

Familiar floating point hardware had provided 36 bits and 60 bits loomed. Experience, thitherto and subsequent, had proven merit to longer floating representations. More floating point precision seemed to be an elixir that solved all sorts of problems that more analysis might have solved in other ways. Longer floating point words were coming anyway and so the analysis was seldom done.

My short foray into small fixed point precision was thus somewhat isolated. (Hydrodynamics had been done with 36 bit fixed point numbers only a few years earlier on machines such as the IBM 701.)

The Physics Model

Today the math would be called “finite element analysis”. Then we called it “difference equations” in contrast to “differential equations” if we called it anything at all. I shall talk of two dimensional Lagrangian methods here where the mesh follows the fluid. Time is explicitly included. The computation is thus three dimensional including time. A mesh, composed of quadrilateral or triangular zones, with rectangular topology covers the region of interest.

An array of small fixed point numbers would hold the x-coordinates of successive vertices in the grid. Another array would hold the y’s. Two more arrays would hold the velocities. Yet other arrays would hold the mass and a few thermodynamical quantities, one per zone. The computation proceeded conceptually by computing the vector separation between vertices of each rectangle, computing the volume, doing some thermodynamics to compute pressure, applying F=ma on each side of the zone thus computing a delta velocity contributed to each vertex by each of the zones which shared the vertex, computing the new vertex velocity, and finally adjusting the coordinates of the vertex according to its velocity. These steps were iterated as computed time and computer time progressed.

This is at least how hydrodynamics had been done in the past. The plan was awkward for the Harvest because the vertex coordinates required too much precision. It seemed that all the physics cared about was the shape of the zone and not its location, except to locate the zone relative to its neighbors.

The answer was easy: store the differences between neighboring vertices. The real vertex positions were needed only when pictures were to be plotted. The real positions were easily computed by summing the differences.

There was a small problem that adding up the delta positions around a path in the mesh should always produce a zero vector. This gave fixed point an advantage because it was easy to arrange the computations to maintain this property mathematically only requiring associativity and commutativity of the arithmetic. (Floating point is not quite associative.)

The same trick was applied to vertex velocities because again the physics does not depend on absolute velocities, only relative velocities. Arrays of relative velocities were kept instead and absolute velocities, when rarely needed could be quickly computed by summation.

At this point the computation had been transformed into a purely Galilean relativistic form. The elements of the computation most demanding of precision had been removed. Only Machian relativity of direction was missing. I did not find a way to capture that.

Much more recently I realized that the lengths of edges of a tetrahedral zone determine its shape and that a computation could proceed carrying only those elements. Also required is a fixed record of the topology of the zones. An author, Tullio Regge, I think, wrote about an n-space being divided into simplexes. Riemann curvature in such a space is concentrated on the n−2 dimensional simplexes. He called them curvature bones. The Bianchi identities provide something like Kerchoff relations for these bones. There is also an extrinsic curvature as it is embedded in 4-space. The extrinsic curvature is a symmetric tensor.