(**************************************************************************************)
(*										                                              *)
(* Copyright 2013 Eloisa Bentivegna						                              *)
(*										                                              *)
(* Analytic.m is a simple Kranc script used to create grid functions which will be    *)
(* used as coefficients (or initial guesses, or exact solutions) by CT_MultiLevel.    *)
(* It is distributed under the General Public License, version 3 or higher.           *)
(*										      *)
(* The runmath.sh and copy-if-changed.sh have been adapted from the same-name         *)
(* scripts in McLachlan; please see McLachlan's licensing and copyright notices       *)
(* regarding this material.                                                           *)
(*										                                              *)
(**************************************************************************************)

Get["KrancThorn`"];


SetEnhancedTimes[False];
(* SetSourceLanguage["C"]; *)

(******************************************************************************)
(* Options *)
(******************************************************************************)

createCode[derivOrder_] :=
Module[{},

prefix = "CT_";
suffix =
  ""
  <> If [derivOrder!=4, "_O" <> ToString[derivOrder], ""]
  ;

ThornName = prefix <> "AnalyticBrillWave" <> suffix;
(* So ThornName = "CT_AnalyticBrillWave" unless I make derivOrder say 6 or 8 later *)

(******************************************************************************)
(* Derivatives *)
(******************************************************************************)

KD = KroneckerDelta;

derivatives =
{
  PDstandardNth[i_]    -> StandardCenteredDifferenceOperator[1,derivOrder/2,i],
  PDstandardNth[i_,i_] -> StandardCenteredDifferenceOperator[2,derivOrder/2,i],
  PDstandardNth[i_,j_] -> StandardCenteredDifferenceOperator[1,derivOrder/2,i] *
                          StandardCenteredDifferenceOperator[1,derivOrder/2,j]
};

FD  = PDstandardNth;


(******************************************************************************)
(* Tensors *)
(******************************************************************************)

(* Register the tensor quantities with the TensorTools package *)
Map [DefineTensor,
     {testinipsi, testinixx, testinixy, testinixz,
      testcxx, testcxy, testcxz, testcyy, testcyz, testczz,
      testcx, testcy, testcz,
      testc0, testc1, testc2, testc3, testc4,
      testXx, testXy, testXz, testAxx, testAxy, testAxz, testAyy, testAyz, testAzz,
      testgxx, testgxy, testgxz, testgyy, testgyz, testgzz, testdetg, testgupxx, 
      testgupxy, testgupyy, testgupzz, x, y, z, gxx, gxy, gxz, gyy, gyz, gzz, detg, gupxx, 
      gupxy, gupxz, gupyy, gupyz, gupzz, Fupx, Fupy, Fupz, g, k, g11, g12, g13, g22, g23, 
      g33, Ricci }];

Map[AssertSymmetricDecreasing,
{
  k[la,lb], g[la,lb]
}];

(******************************************************************************)
(* Shorthands *)
(******************************************************************************)



(* Use the CartGrid3D variable names *)
x1=x; x2=y; x3=z;

(******************************************************************************)
(* Expressions *)
(******************************************************************************)

detgExpr  -> Det [MatrixOfComponents [g [la,lb]]];
one -> 1;

(******************************************************************************)
(* Groups *)
(******************************************************************************)

evolvedGroups = {};
evaluatedGroups =
  {SetGroupName [CreateGroupFromTensor [testinipsi], prefix <> "testinipsi"],
   SetGroupName [CreateGroupFromTensor [testcxx  ], prefix <> "testcxx"],
   SetGroupName [CreateGroupFromTensor [testcxy  ], prefix <> "testcxy"],
   SetGroupName [CreateGroupFromTensor [testcxz  ], prefix <> "testcxz"],
   SetGroupName [CreateGroupFromTensor [testcyy  ], prefix <> "testcyy"],
   SetGroupName [CreateGroupFromTensor [testcyz  ], prefix <> "testcyz"],
   SetGroupName [CreateGroupFromTensor [testczz  ], prefix <> "testczz"],
   SetGroupName [CreateGroupFromTensor [testcx   ], prefix <> "testcx"],
   SetGroupName [CreateGroupFromTensor [testcy   ], prefix <> "testcy"],
   SetGroupName [CreateGroupFromTensor [testcz   ], prefix <> "testcz"],
   SetGroupName [CreateGroupFromTensor [testc0   ], prefix <> "testc0"],
   SetGroupName [CreateGroupFromTensor [testc1   ], prefix <> "testc1"],
   SetGroupName [CreateGroupFromTensor [testc2   ], prefix <> "testc2"],
   SetGroupName [CreateGroupFromTensor [testc3   ], prefix <> "testc3"],
   SetGroupName [CreateGroupFromTensor [testc4   ], prefix <> "testc4"],
   SetGroupName [CreateGroupFromTensor [testXx   ], prefix <> "testXx"],
   SetGroupName [CreateGroupFromTensor [testXy   ], prefix <> "testXy"],
   SetGroupName [CreateGroupFromTensor [testXz   ], prefix <> "testXz"],
   SetGroupName [CreateGroupFromTensor [testAxx  ], prefix <> "testAxx"],
   SetGroupName [CreateGroupFromTensor [testAxy  ], prefix <> "testAxy"],
   SetGroupName [CreateGroupFromTensor [testAxz  ], prefix <> "testAxz"],
   SetGroupName [CreateGroupFromTensor [testAyy  ], prefix <> "testAyy"],
   SetGroupName [CreateGroupFromTensor [testAyz  ], prefix <> "testAyz"],
   SetGroupName [CreateGroupFromTensor [testAzz  ], prefix <> "testAzz"],
   SetGroupName [CreateGroupFromTensor [testgxx  ], prefix <> "testgxx"],
   SetGroupName [CreateGroupFromTensor [testgxy  ], prefix <> "testgxy"],
   SetGroupName [CreateGroupFromTensor [testgyy  ], prefix <> "testgyy"],
   SetGroupName [CreateGroupFromTensor [testgzz  ], prefix <> "testgzz"],
   SetGroupName [CreateGroupFromTensor [testdetg  ], prefix <> "testdetg"],
   SetGroupName [CreateGroupFromTensor [testgupxx  ], prefix <> "testgupxx"],
   SetGroupName [CreateGroupFromTensor [testgupxy  ], prefix <> "testgupxy"],
   SetGroupName [CreateGroupFromTensor [testgupyy  ], prefix <> "testgupyy"],
   SetGroupName [CreateGroupFromTensor [testgupzz  ], prefix <> "testgupzz"]
  };

declaredGroups = Join [evolvedGroups, evaluatedGroups];

declaredGroupNames = Map [First, declaredGroups];



extraGroups =
  {{"Grid::coordinates", {x, y, z}}
};

admGroups =
  {{"admbase::metric", {gxx,gxy,gxz,gyy,gyz,gzz}},
   {"admbase::curv", {kxx,kxy,kxz,kyy,kyz,kzz}},
   {"admbase::lapse", {alp}},
   {"admbase::shift", {betax,betay,betaz}}};



groups = Join [declaredGroups, extraGroups, admGroups];

(******************************************************************************)
(* Metric and Extrinsic Curvature Components *)
(******************************************************************************)

Axx[x_, y_, z_] := 0;
Axy[x_, y_, z_] := 0;
Axz[x_, y_, z_] := 0;
Ayy[x_, y_, z_] := 0;
Ayz[x_, y_, z_] := 0;
Azz[x_, y_, z_] := 0;

epsi -> one,
rho1 -> Sqrt[x*x + y*y],
e2q -> Exp[2*q[rho1,z]],
rho2 -> x*x + y*y,
gxx -> e2q + (one - e2q)*y*y/rho2,
gxy -> - (one - e2q)*x*y/rho2,
gyy -> e2q + (one - e2q)*x*x/rho2,
gzz -> e2q,
detg -> gxx gyy gzz - gxy gxy gzz,
oneoverdetg -> 1/detg,
gupxx -> 1/gxx + gxy gxy/(gxx gxx gyy -gxx gxy gxy),
gupxy -> -gxy/(gxx gyy -gxy gxy),
gupyy -> 1/(gyy - gxy gxy / gxx),
gupzz -> 1/gzz,

g11 = gxx;
g12 = gxy;
g13 = gxz;
g22 = gyy;
g23 = gyz;
g33 = gzz;


(* Brill q function *)
q[rho1_,zcyl_] := -ampBrill*rho1*((rho1/sigmarho)^2 + (zcyl/sizmaz)^2),

TrRicciCalc =
{
  Name -> TrRicci,
  Schedule -> {"AT CCTK_INITIAL before CT_MultiLevel"},
  Where -> Interior,
  Shorthands -> {detg, gu[ua,ub], G[ua,lb,lc], Ricci[la,lb]},
  Equations ->
  {
    detg -> detgExpr,
    gu[ua,ub] -> 1/detg detgExpr MatrixInverse[g[ua,ub]], 
    G[ua,lb,lc] -> 1/2 gu[ua,ud]
                   (PD[g[lb,ld],lc] + PD[g[lc,ld],lb] - PD[g[lb,lc],ld]),
    Ricci[la,lb] -> gu[us,ur](G[um,la,lr] G[uk,ls,lb] g[lk,lm] - G[um,la,lb] G[uk,ls,lr] g[lk,lm])
                + 1/2 gu[uc,ud] (- PD[g[lc,ld],la,lb] + PD[g[lc,la],ld,lb]
                                 - PD[g[la,lb],lc,ld] + PD[g[ld,lb],lc,la]),
    trRicci -> Ricci[la,lb] gu[ua,ub] (*  trRicci to be used in forming coeff c0 = -1/8 trRicci below *)
  }
};


AnalyticBrillWaveCalc =
{
  Name -> ThornName ,
  Schedule -> {"AT CCTK_INITIAL before CT_MultiLevel"},
(*  ConditionalOnKeyword -> {"free_data", "BrillWave_K"}, *)
  Where -> Everywhere,
  Equations ->
  {
        sqrtdetg -> Sqrt[detg],
        oneoversqrtdetg -> 1/sqrtdetg,
	Fupx -> oneoversqrtdetg ( D[sqrtdetg gupxx,x] + D[sqrtdetg gupxy,y] +  D[sqrtdetg gupxz,z]),
	Fupy -> oneoversqrtdetg ( D[sqrtdetg gupxy,x] + D[sqrtdetg gupyy,y]  + D[sqrtdetg gupxz,z]),
	Fupz -> oneoversqrtdetg ( D[sqrtdetg gupxz,x] + D[sqrtdetg gupyz,y]  + D[sqrtdetg gupzz,z]),
	coeffpartialwrtxx = gupxx,
	coeffpartialwrtxy = gupxy,
        coeffpartialwrtyy = gupyy,
	coeffpartialwrtzz = gupzz,
	coeffpartialwrtx  = Fupx,
	coeffpartialwrty  = Fupy,
        coeffpartialwrtz  = Fupz,
    testinipsi -> 1.0, 
    testcxx ->  coeffpartialwrtxx, 
    testcxy ->  coeffpartialwrtxy, 
    testcyy ->  coeffpartialwrtyy, 
    testczz ->  coeffpartialwrtzz, 
    testcx  ->   coeffpartialwrtx,
    testcy  ->   coeffpartialwrty,
    testcz  ->   coeffpartialwrtz,
    testc0  ->  -(1/8)*trRicci,
    testXx  -> 0.0,
    testXy  -> 0.0,
    testXz  -> 0.0,
    testAxx -> 0.0,
    testAxy -> 0.0,
    testAxz -> 0.0,
    testAyy -> 0.0,
    testAyz -> 0.0,
    testAzz -> 0.0
   }
};





(******************************************************************************)
(* Implementations *)
(******************************************************************************)

(******************************************************************************)
(* Parameters *)
(******************************************************************************)

intParameters =
{
}

realParameters =
{

  {
    Name -> ampBrill,
    Description -> "Coefficient in the q function in initial 3-metric",
    Default -> 0
  },

  {  Name -> sigmarho,
    Description -> " rho width of Gaussian in q function in initial 3-metric",
    Default -> 1
  },
  
  {  Name -> sigmaz,
    Description -> " z width of Gaussian in q function in initial 3-metric",
    Default -> 1
  },
  {
    Name -> eps,
    Description -> "Smoothing factor",
    Default -> 1`*^-6
  }
 
};

keywordParameters =
{
  {
    Name -> "free_data",
    Description -> "How to set the free data for the extrinsic curvature?",
    AllowedValues -> {"BrillWave_K"},
    Default -> "BrillWave_K"
  }
};

(******************************************************************************)
(* Construct the thorns *)
(******************************************************************************)

calculations =
{
(*  TrRicciCalc, *)
  AnalyticBrillWaveCalc
};

CreateKrancThornTT [groups, ".", ThornName,
  Calculations -> calculations,
  DeclaredGroups -> declaredGroupNames,
  PartialDerivatives -> derivatives,
  UseLoopControl -> True,
  UseVectors -> False,
  IntParameters -> intParameters,
  RealParameters -> realParameters,
  InheritedImplementations -> {"admbase"},
  KeywordParameters -> keywordParameters
];

]


(******************************************************************************)
(* Options *)
(******************************************************************************)

(* derivative order: 2, 4, 6, 8, ... *)
(* useGlobalDerivs: False or True *)
(* timelevels: 2 or 3
   (keep this at 3; this is better chosen with a run-time parameter) *)
(* matter: 0 or 1
   (matter seems cheap; it should be always enabled) *)

createCode[4];

