Compute covariant derivative in Mathematica

Here is some code I wrote a while back:

ClearAll["Global`*"]
SetAttributes[Rs, Constant]
$Assumptions = Rs > 0;
Coordinates = {t, r, \[Theta], \[Phi]};
dim = Length[Coordinates];

MetricTensorLL = {{(1 - Rs/r), 0, 0, 0}, {0, -(1 - Rs/r)^-1, 0, 
 0}, {0, 0, -r^2, 0}, {0, 0, 0, -r^2 Sin[\[Theta]]^2}};

MetricTensorUU := MetricTensorUU = Simplify[Inverse[MetricTensorLL]];

ChristoffelSymbolsULL := 
 ChristoffelSymbolsULL = 
  Simplify@
   Array[1/2 Sum[
    MetricTensorUU[[#1, \[Lambda]]] (D[
     MetricTensorLL[[#3, \[Lambda]]], Coordinates[[#2]]] + 
      D[MetricTensorLL[[\[Lambda], #2]], Coordinates[[#3]]] - 
       D[MetricTensorLL[[#2, #3]], 
        Coordinates[[\[Lambda]]]]), {\[Lambda], dim}] &, {dim, 
         dim, dim}];

ChristoffelSymbolsLLL := 
 ChristoffelSymbolsLLL = 
  Simplify@
   Array[1/2 (D[MetricTensorLL[[#1, #2]], Coordinates[[#3]]] + 
    D[MetricTensorLL[[#1, #3]], Coordinates[[#2]]] - 
     D[MetricTensorLL[[#2, #3]], Coordinates[[#1]]]) &, {dim, dim,
      dim}];

RiemannCurvatureTensorULLL := 
 RiemannCurvatureTensorULLL = 
  Simplify@
   Array[D[ChristoffelSymbolsULL[[#1, #2, #4]], Coordinates[[#3]]] - 
    D[ChristoffelSymbolsULL[[#1, #2, #3]], Coordinates[[#4]]] + 
     Sum[ChristoffelSymbolsULL[[#1, #3, \[Epsilon]]] \
      ChristoffelSymbolsULL[[\[Epsilon], #2, #4]], {\[Epsilon], dim}] - 
       Sum[ChristoffelSymbolsULL[[#1, #4, \[Epsilon]]] \
        ChristoffelSymbolsULL[[\[Epsilon], #2, #3]], {\[Epsilon], 
         dim}] &, {dim, dim, dim, dim}];

RiemannCurvatureTensorLLLL := 
 RiemannCurvatureTensorLLLL = 
  Simplify@
   Array[Sum[
    MetricTensorLL[[#1, \[Tau]]] \
     RiemannCurvatureTensorULLL[[\[Tau], #2, #3, #4]], {\[Tau], 
      dim}] &, {dim, dim, dim, dim}];

RiemannCurvatureTensorUUUU := 
 RiemannCurvatureTensorUUUU = 
  Simplify@
   Array[Sum[
    MetricTensorUU[[#2, \[Alpha]]] MetricTensorUU[[#3, \[Beta]]] \
     MetricTensorUU[[#4, \[Gamma]]] RiemannCurvatureTensorULLL[[#1, \
      \[Alpha], \[Beta], \[Gamma]]], {\[Alpha], dim}, {\[Beta], 
       dim}, {\[Gamma], dim}] &, {dim, dim, dim, dim}];

RiemannCurvatureTensorLL := 
 RiemannCurvatureTensorLL = 
  Simplify@
   Array[Sum[
    RiemannCurvatureTensorULLL[[\[Lambda], #1, \[Lambda], #2]], {\
     \[Lambda], dim}] &, {dim, dim}];

RiemannCurvatureTensorUL := 
 RiemannCurvatureTensorUL = 
  Simplify@
   Array[Sum[
    MetricTensorUU[[#1, \[Lambda]]] RiemannCurvatureTensorLL[[\
     \[Lambda], #2]], {\[Lambda], dim}] &, {dim, dim}];

ScalarCurvature := ScalarCurvature = Tr[RiemannCurvatureTensorUL];

KretschmannScalar := 
 KretschmannScalar = 
  Simplify@
   Sum[RiemannCurvatureTensorLLLL[[\[Alpha], \[Beta], \[Gamma], \
    \[Delta]]] RiemannCurvatureTensorUUUU[[\[Alpha], \[Beta], \[Gamma], \
     \[Delta]]], {\[Alpha], dim}, {\[Beta], dim}, {\[Gamma], 
      dim}, {\[Delta], dim}];

WeylCurvatureTensorLLLL := 
 WeylCurvatureTensorLLLL = 
  Simplify@
   Array[If[dim > 3, 
    RiemannCurvatureTensorLLLL[[#1, #2, #3, #4]] - 
     1/(
      dim - 2) (MetricTensorLL[[#1, #3]] \
       RiemannCurvatureTensorLL[[#4, #2]] + 
        MetricTensorLL[[#2, #4]] RiemannCurvatureTensorLL[[#3, \
         #1]] - MetricTensorLL[[#1, #4]] RiemannCurvatureTensorLL[[#3, #2]] - MetricTensorLL[[#2, #3]] RiemannCurvatureTensorLL[[#4, \
          #1]]) + ScalarCurvature/((dim - 1) (dim - 2)) (MetricTensorLL[[#1, #3]] MetricTensorLL[[#4, #2]] - 
           MetricTensorLL[[#1, #4]] MetricTensorLL[[#3, #2]]), 
            0] &, {dim, dim, dim, dim}];

EinsteinTensor := 
 EinsteinTensor = 
  Simplify[
   RiemannCurvatureTensorLL - 1/2 MetricTensorLL ScalarCurvature];

ConformallyFlatSpaceQ := 
 ConformallyFlatSpaceQ = 
  Simplify[Equal[Sequence @@ Flatten@WeylCurvatureTensorLLLL, 0]];

MaximallySymmetricSpaceQ := 
 MaximallySymmetricSpaceQ = 
  Simplify[
   And @@ Flatten@
    Map[# == 0 &, 
     RiemannCurvatureTensorLLLL - 
      Array[ScalarCurvature/(
       dim (dim - 
        1)) (MetricTensorLL[[#1, #3]] MetricTensorLL[[#2, #4]] -
          MetricTensorLL[[#1, #4]] MetricTensorLL[[#2, #3]]) &, \
           {dim, dim, dim, dim}], {4}]];

It is hopefully self-explanatory.

In this particular example I am using the Schwarzschild metric. The first line declares Rs, the Schwarzschild radius, to be a constant. We then define the coordinates of the manifold, in this case, $t,r,\theta,\phi$. After that, we input the components of the metric tensor, with its indices lowered (as indicated by the LL tag, meaning "lower-lower").

The rest of tensors are calculated by MMA. For example, ChristoffelSymbolsULL are the components of the Christoffel symbols, with one upper and two lower indices. Similarly, ChristoffelSymbolsLLL are the components of the Christoffel symbols with all indices lowered.

The code also computes the Riemann tensor, the Weyl tensor, and several scalars (Ricci and Kretschmann). Finally, there is a check for whether the manifold is conformally flat and/or maximally symmetric.

To compute covariant derivatives, you can use the known value of the Christoffel symbols, or the expression

Sum[1/Sqrt[Det[MetricTensorLL]] D[Sqrt[Det[MetricTensorLL]] MetricTensorUU[[\[Mu] + 1, \[Nu] + 1]] D[f @@ Coordinates, Coordinates[[\[Nu] + 1]]], Coordinates[[\[Mu] + 1]]], {\[Mu], 0, dim - 1}, {\[Nu], 0, dim - 1}]

where f is an arbitrary scalar function. For higher rank tensor fields, you'll have to make some small modifications to the code.

There is undocumented functionality in the SymbolicTensors package, which underlies CoordinateChartData, CoordinateTransformData, and the coordinate-system awareness of Grad, etc. You could add the context to $ContextPath save needing to type everywhere, though I won't do that in the example below.

First, you want to define a "patch". This can be done in various ways, but for a diagonal metric tensor the easiest way is scale factors. Here I'll use Minkowski in spherical coordinates as my example:

vars = {t, r, \[Theta], \[Phi]};
patch = SymbolicTensors`ScaleFactorGeometryPatch[{-1, 1, r, r Sin[\[Theta]]}, vars];

This patch is the equivalent of CoordinateChartData for your custom metric. Take a look at patch["Properties"] for some interesting things it will spit out.

Next, we need to define the two vector fields in the tensor language of the package. I'm just using the SlotSequence notation to save myself some typing, The evaluated form would work just as well for input. Below, TangentBasis is a representation of the coordinated basis for vector fields. There is a CotangentBasis for the space of one-forms. If you wish to use, e.g., an orthonormal basis, the TangentBasis with TransformedBasis and supply the change of basis matrix in the second argument.

v = SymbolicTensors`Tensor[{vt[##], vr[##], v\[Theta][##], 
     v\[CurlyPhi][##]}, {SymbolicTensors`TangentBasis[{##}]}] & @@ vars

u = SymbolicTensors`Tensor[{ut[##], ur[##], u\[Theta][##], 
     u\[CurlyPhi][##]}, {SymbolicTensors`TangentBasis[{##}]}] & @@ vars

Finally, we want to do an actual computation. For this there is a CovariantD which returns $\nabla_b v^a$, using the abstract index notation. CovariantD has the same sytnax as Grad and friends, and all of these take a patch in the third argument. (Note that Grad returns the raised derivative $\nabla^bv^a$.) The covariant derivative along a vector fields is simply $u^b\nabla_bv^a$. Since all differentiation functions in M- add the new slot at the end, this then is simply:

SymbolicTensors`CovariantD[v, vars, patch] . u

If you don't like the fact u comes at the end, you could write this more verbosely as

TensorContract[ 
    TensorProduct[u, SymbolicTensors`CovariantD[v, vars, patch]],
   {1, 3}
]