Improve xTensor support to the point where it can generate an xTensor version of the EM example. There are still several very hack parts which will need to be improved.

2 files changed, 177 insertions, 21 deletions

diff --git a/Examples/EM-xTensor.m b/Examples/EM-xTensor.m
new file mode 100644
index 0000000..944bed3
--- /dev/null
+++ b/Examples/EM-xTensor.m
@@ -0,0 +1,156 @@
+(* ::Package:: *)
+
+(*  Copyright 2004 Sascha Husa, Ian Hinder, Christiane Lechner, Barry Wardell
+
+    This file is part of Kranc.
+
+    Kranc is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    Kranc is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Foobar; if not, write to the Free Software
+    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+*)
+
+$KrancTensorPackage="xTensorKranc`";
+Get["KrancTensor`"];
+Get["GeneralRelativity`LoadMetric`"];
+$DefInfoQ=False;
+$CVVerbose=False;
+$PrePrint=ScreenDollarIndices;
+
+(**************************************************************************************)
+(* Derivatives *)
+(**************************************************************************************)
+
+derivatives =
+{
+  PDstandard2nd[i_] -> StandardCenteredDifferenceOperator[1,1,i],
+  PDstandard2nd[i_, i_] -> StandardCenteredDifferenceOperator[2,1,i],
+  PDstandard2nd[i_, j_] -> StandardCenteredDifferenceOperator[1,1,i] *
+    StandardCenteredDifferenceOperator[1,1,j],
+
+  PDstandard4th[i_] -> StandardCenteredDifferenceOperator[1,2,i],
+  PDstandard4th[i_, i_] -> StandardCenteredDifferenceOperator[2,2,i],
+  PDstandard4th[i_, j_] -> StandardCenteredDifferenceOperator[1,2,i] *
+    StandardCenteredDifferenceOperator[1,2,j]
+};
+
+(**************************************************************************************)
+(* Tensors *)
+(**************************************************************************************)
+(* Load Euclidean metric and set components of g, epsilon and gamma *)
+LoadMetric["Euclidean"];
+MetricCompute[g, Euclidean, "Christoffel"[1,-1,-1], CVSimplify -> Simplify];
+AllComponentValues[ToBasis[Euclidean][epsilong[-a,e,f]],
+  Table[Signature[{i,j,k}],{i,3},{j,3},{k,3}]];
+RuleToSet[ChristoffelCDPDEuclidean];
+RuleToSet[g];
+RuleToSet[epsilong];
+
+(* Register the tensor quantities with the TensorTools package *)
+DefTensor[El[-a],EuclideanM];
+DefTensor[B[-a],EuclideanM];
+
+(**************************************************************************************)
+(* Groups *)
+(**************************************************************************************)
+
+(* NB This will give B the symmetries of a VECTOR, which is not correct *)
+evolvedGroups = Map[CreateGroupFromTensor, ToBasis[Euclidean][{El[-a], B[-a]}]];
+evaluatedGroups = 
+ {{"constraints", {CEl, CB}},
+  {"endens",{rho}}};
+
+declaredGroups = Join[evolvedGroups, evaluatedGroups];
+declaredGroupNames = Map[First, declaredGroups];
+
+groups = Join[declaredGroups];
+
+(**************************************************************************************)
+(* Initial data *)
+(**************************************************************************************)
+
+initialCalc =
+{
+  Name -> "EM_initial",
+  Schedule -> {"at CCTK_INITIAL"},
+  Equations -> 
+  {
+    El1 -> sigma*Cos[2 Pi (x + y)] ,
+    El2 -> - (1 - sigma)*Cos[2 Pi x] - sigma*Cos[2 Pi (x + y)],
+    El3 -> 0,
+    B1  -> 0,
+    B2  -> 0,
+    B3  -> (1 - sigma)*Cos[2 Pi x] + sigma*Cos[2 Pi (x + y)]
+  }
+};
+
+(**************************************************************************************)
+(* Evolution equations *)
+(**************************************************************************************)
+
+evolCalc =
+{
+  Name -> "EM_evol",
+  Schedule -> {"in MoL_CalcRHS"},
+  Equations -> 
+  {
+    dot[ToBasis[Euclidean][El[-a]]] -> ToBasis[Euclidean][(epsilong[-a,e,f]  CD[-e][B[-f]])],
+    dot[ToBasis[Euclidean][ B[-a]]] -> ToBasis[Euclidean][-(epsilong[-a,e,f] CD[-e][El[-f]])]
+  }
+};
+
+(**************************************************************************************)
+(* Constraint equations *)
+(**************************************************************************************)
+
+constraintsCalc =
+{
+  Name -> "EM_constraints",
+  Equations -> 
+  {
+    CEl -> ToBasis[Euclidean][g[a,b]CD[-b][El[-a]]],
+    CB -> ToBasis[Euclidean][ g[a,b]CD[-b][B[-a]]] 
+  }
+};
+
+(**************************************************************************************)
+(* Energy equation *)
+(**************************************************************************************)
+
+energyCalc =
+{
+  Name -> "EM_energy",
+  Equations -> 
+  {
+    rho -> ToBasis[Euclidean][El[-a] El[a]/2 + B[-a] B[a]/2]
+  }
+};
+
+realParameters = {sigma};
+
+(**************************************************************************************)
+(* Construct the thorn *)
+(**************************************************************************************)
+
+calculations = 
+{
+  initialCalc,
+  evolCalc, 
+  constraintsCalc, 
+  energyCalc
+};
+
+CreateKrancThornTT[groups, ".", "EM-xTensor", 
+  Calculations -> calculations,
+  DeclaredGroups -> declaredGroupNames,
+  PartialDerivatives -> derivatives,
+  RealParameters -> realParameters];
diff --git a/Tools/CodeGen/xTensorKranc.m b/Tools/CodeGen/xTensorKranc.m
index 2ee92d6..f92072f 100644
--- a/Tools/CodeGen/xTensorKranc.m
+++ b/Tools/CodeGen/xTensorKranc.m
@@ -17,7 +17,7 @@
     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
 *)
 
-BeginPackage["xTensorKranc`", {"Kranc`", "xAct`xTensor`", "xAct`xCore`", "xAct`xCoba`"}];
+BeginPackage["xTensorKranc`", {"Differencing`", "Kranc`", "KrancGroups`", "xAct`xTensor`", "xAct`xCore`", "xAct`xCoba`"}];
 
 CreateGroupFromTensor::usage = "";
 ReflectionSymmetries::usage = "Produce a list of reflection symmetries of a tensor.";
@@ -34,35 +34,34 @@ TensorCharacterString[t_Symbol?xTensorQ[inds___]] := StringJoin[If[UpIndexQ[#],"
 Options[makeExplicit] = {IncludeCharacter -> False};
 SetAttributes[makeExplicit, Listable];
 e : makeExplicit[_Plus, opts:OptionsPattern[]] := Distribute[Unevaluated[e]];
-makeExplicit[x_Times, opts:OptionsPattern[]] :=Map[makeExplicit[#, opts]&, x];
+makeExplicit[x_Times, opts:OptionsPattern[]] := Map[makeExplicit[#, opts]&, x];
 makeExplicit[Power[x_,p_], opts:OptionsPattern[]] := Power[makeExplicit[x, opts],p];
-makeExplicit[x_?NumericQ, OptionsPattern[]] := x;
-makeExplicit[x_, {}] := makeExplicit[x];
-makeExplicit[t_Symbol?xTensorQ[inds___], OptionsPattern[]] := Module[{indexNumbers,character,indexString},
+makeExplicit[x_?NumericQ, opts:OptionsPattern[]] := x;
+makeExplicit[x_, {}, opts:OptionsPattern[]] := makeExplicit[x];
+makeExplicit[t_Symbol?xTensorQ[inds___], opts:OptionsPattern[]] := Module[{indexNumbers,character,indexString},
   indexNumbers=First/@{inds};
   
   If[OptionValue[IncludeCharacter],
     character = TensorCharacterString[t[inds]];
     indexString = StringJoin[character,ToString/@indexNumbers],
-
     indexString = StringJoin[ToString/@indexNumbers]
   ];
   SymbolJoin[PrintAs[t],Sequence@@indexString]
 ];
 
-makeExplicit[cd_?CovDQ[ind_][expr_], OptionsPattern[]] := Module[{indexNumbers,character,indexString},
-  indexNumbers=First/@{ind};
+makeExplicit[x_, opts:OptionsPattern[]] := x;
 
-  If[OptionValue[IncludeCharacter],
-    character = TensorCharacterString[cd[ind]];
-    indexString = StringJoin[character,ToString/@indexNumbers],
+makeExplicit[(cd_?CovDQ)[ind_][expr_], opts:OptionsPattern[]] := Module[{indexNumbers},
+  indexNumbers=First/@{ind};
 
-    indexString = StringJoin[ToString/@indexNumbers]
-  ];
-  SymbolJoin["d",indexString,makeExplicit[expr]]
+  Global`PDstandard2nd[makeExplicit[expr, opts], Sequence@@indexNumbers]
 ];
 
 Options[ExpandComponents] = Options[ExpandComponents];
+
+ExpandComponents[x_Rule, opts:OptionsPattern[makeExplicit]] := Thread[ExpandComponents[x[[1]], opts] -> ExpandComponents[x[[2]], opts]];
+ExpandComponents[dot[x_], opts:OptionsPattern[makeExplicit]] := dot/@ExpandComponents[x, opts];
+ExpandComponents[x_List, opts:OptionsPattern[makeExplicit]] := Flatten[Map[ExpandComponents[#, opts]&, x], 1];
 ExpandComponents[x_, opts:OptionsPattern[makeExplicit]] :=
   Module[{eqs, options},
 
@@ -74,9 +73,10 @@ ExpandComponents[x_, opts:OptionsPattern[makeExplicit]] :=
   ]
 ];
 
+
 (* Compute the reflection symmetries of a tensor *)
-ReflectionSymmetries[t_Symbol?xTensorQ[inds__], b_] :=
-  Module[{cnums, components, componentIndices, counts},
+ReflectionSymmetries[t_Symbol?xTensorQ[inds__]] :=
+  Module[{b=Global`Euclidean, cnums, components, componentIndices, counts},
     (* Get the compoent indices of the basis *)
     cnums = CNumbersOf[b, VBundleOfBasis[b]];
 
@@ -86,19 +86,19 @@ ReflectionSymmetries[t_Symbol?xTensorQ[inds__], b_] :=
     (* Get the indices of each component *)
     componentIndices = Map[IndicesOf[b], components];
 
-
     (* Count the number of instances of each basis index. *)
     countInds[expr_, basis_, cinds_] := Map[(Count[expr,{#,basis}]+Count[expr,{#,-basis}])&, cinds];
     counts = Map[countInds[#, b, cnums]&, componentIndices];
 
     (* For each instance, multiply by -1 *)
-    Thread[components -> (-1)^counts]
+    Thread[ExpandComponents[t[inds]] -> (-1)^counts]
 ];
 
-ReflectionSymmetries[t_Symbol?xTensorQ[], b_] := {};
+ReflectionSymmetries[t_Symbol?xTensorQ[]] := t -> {1,1,1};
+ReflectionSymmetries[t_] := t -> {1, 1, 1};
 
 (* FIXME: Implement this fully *)
-GetTensorAttribute[t_Symbol?xTensorQ[inds___], TensorWeight] := WeightOfTensor[t];
+GetTensorAttribute[t_Symbol?xTensorQ, TensorWeight] := WeightOfTensor[t];
 
 CreateGroupFromTensor[t_Symbol?xTensorQ[inds__]] := Module[{tCharString, nInds, tags, vars, group},
   InfoMessage[InfoFull, "Creating group from tensor with kernel " <> SymbolName[t] <> " and indices " <> ToString[{inds}]];
@@ -120,7 +120,7 @@ CreateGroupFromTensor[t_Symbol?xTensorQ[inds__]] := Module[{tCharString, nInds,
   Return[group]
 ];
 
-ReflectionSymmetries[x___]:= Throw["ReflectionSymmetriesOfTensor error: "<>ToString[x]];
+ReflectionSymmetries[x___]:= Throw["ReflectionSymmetries error: "<>ToString[x]];
 CreateGroupFromTensor[x___]:= Throw["CreateGroupFromTensor error: "<>ToString[x]];
 
 CheckTensors[expr_] := Validate[expr];