(* ::Package:: *)

(* ::Input:: *)
(*$Path=Join[$Path,{"../../../repos/Kranc/Tools/CodeGen","../../../repos/Kranc/Tools/MathematicaMisc"}];*)
(**)
(*Get["KrancThorn`"];*)
(**)
(*(*SetDebugLevel[InfoFull];*)*)
(**)
(*SetEnhancedTimes[False];*)
(**)
(*(****************************************************************************Derivatives****************************************************************************)*)
(**)
(*derivOrder=2;*)
(**)
(*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;*)
(*ResetJacobians;*)
(*DefineJacobian[PD,FD,KD,Zero3];*)
(*(****************************************************************************Tensors****************************************************************************)(*Register all the tensors that will be used with TensorTools*)*)
(*Map[DefineTensor,{wavePsi,wavePi}];*)
(*(*Determinants of the metrics in terms of their components (Mathematica symbolic expressions)*)gDet=Det[MatrixOfComponents[g[la,lb]]];*)
(*(******************************************************************************)*)
(*(* Declare tensors *)*)
(*(******************************************************************************)*)
(**)
(*(* Rename some Cactus grid functions *)*)
(**)
(*(* CartGrid3D coordinates *)*)
(*x1=x; x2=y; x3=z;*)
(*g11=gxx; g12=gxy; g22=gyy; g13=gxz; g23=gyz; g33=gzz;*)
(*admalpha=alp;*)
(*beta1=betax; beta2=betay; beta3=betaz;*)
(*(* TmunuBase variables *)*)
(*T00=eTtt;*)
(*T01=eTtx; T02=eTty; T03=eTtz;*)
(*T11=eTxx; T12=eTxy; T22=eTyy; T13=eTxz; T23=eTyz; T33=eTzz;*)
(*DefineTensor[alpha];*)
(*DefineTensor[beta, {ua}];*)
(**)
(*DefineTensor[T00];*)
(*DefineTensor[T0, {la}];*)
(**)
(**)
(**)
(*(***************************************************************************Groups****************************************************************************)*)
(*SetGroupTimelevels[g_,tl_]=Join[g,{Timelevels->tl}];*)
(*evolvedGroups={CreateGroupFromTensor[wavePsi],(*psi*)CreateGroupFromTensor[wavePi](*pi*)};*)
(*evaluatedGroups={(*Tab*)CreateGroupFromTensor[T[la,lb]]};*)
(*declaredGroups=Join[evolvedGroups,evaluatedGroups];*)
(*extraGroups =*)
(*  {{"grid::coordinates", {x, y, z, r}},*)
(*   {"ADMBase::metric",  {gxx, gxy, gxz, gyy, gyz, gzz}},*)
(*   {"ADMBase::lapse",   {alp}},*)
(*   {"ADMBase::shift",   {betax, betay, betaz}},*)
(*   {"TmunuBase::stress_energy_scalar", {eTtt}},*)
(*   {"TmunuBase::stress_energy_vector", {eTtx, eTty, eTtz}},*)
(*   {"TmunuBase::stress_energy_tensor", {eTxx, eTxy, eTxz, eTyy, eTyz, eTzz}}};*)
(*groups = Join[declaredGroups, extraGroups];*)
(*declaredGroupNames=Map[First,groups];*)
(**)
(*groups=Join[declaredGroups];*)
(*(******************************************************************************)*)
(*(*Initial data*)*)
(*(******************************************************************************)*)
(*waveCalc={Name->"waveKG",Schedule->{"IN ADMBase_InitialData"},Equations->{wavePsi->A Exp[-(x^2+y^2+z^2) / sigma^2],wavePi->2 A x vExp[-(x^2+y^2+z^2) / sigma^2] / sigma^2}};*)
(*(******************************************************************************)*)
(*(*Evolution equations*)*)
(*(******************************************************************************)*)
(*evolCalc={*)
(*Name->"waveKG_RHS",*)
(*Schedule->{"IN waveKG_evolCalcGroup"},*)
(*Where->Interior,*)
(*Equations->{*)
(*(*dt psi=pi*)*)
(*dot[wavePsi]->wavePi,*)
(*(*dt pi=2 beta[ua] da pi + alpha^2 g[ua,ub] da db psi - alpha^2 psi*)*)
(*dot[wavePi]->2 beta[ua] PD[wavePi,la] + alpha^2 g[ua,ub] PD[wavePsi, la,lb] - alpha^2 wavePsi,*)
(*(*Assign Tab=Da psi Db psi - 1/2 gab (D[lc]psi D[uc]psi + psi^2)*)*)
(*T[la,lb]-> CDtoPD[CD[wavePsi,la]] CDtoPD[CD[wavePsi,lb]] - (1/2) g[la,lb] (CDtoPD[wavePsi,lc] CDtoPD[wavePsi,uc] + wavePsi^2)}};*)
(**)
(*(******************************************************************************)*)
(*(*Parameters*)*)
(*(******************************************************************************)*)
(**)
(*realParameters={{Name->A,Description->"wave amplitude",Default->0.001},{Name->sigma,Description->"width",Default->1}};*)
(*(***************************************************************************************)*)
(*(*Thorn construction********************************************************************)*)
(*(***************************************************************************************)*)
(*calculations =*)
(*{*)
(*  waveCalc,*)
(*  evolCalc*)
(*};*)
(**)
(*CreateKrancThornTT [groups, ".", "ML_waveKG",*)
(*  Calculations -> calculations,*)
(*  DeclaredGroups -> declaredGroupNames,*)
(*  PartialDerivatives -> derivatives,*)
(*  UseLoopControl -> True,*)
(*  RealParameters -> realParameters*)
(*];*)
(**)
