Kranc source file demonstrating duplicated reads/writes entries in schedule.ccl

CT_BrillAnalytic.m

Clear["`*"];
(* ::Package:: *)

(**************************************************************************************)
(*										                                              *)
(* 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`"];

Print["Just after Get KrancThorn "];
SetEnhancedTimes[False];
(* SetSourceLanguage["C"]; *)

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

createCode[derivOrder_] :=
Module[{},

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

ThornName = prefix <> "BrillAnalytic" <> suffix;
(* So ThornName = "CT_BrillAnalytic" 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]
};

PD = 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,
      x, y, z, detg,
      gf, FF, g, k, Ricci, G, trRicci, gu}];

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

Map [AssertSymmetricIncreasing,
     {g[la,lb], k[la,lb]}];
Map [AssertSymmetricIncreasing,
     {gu[ua,ub]}];
	 (* CD added next *)
Map [AssertSymmetricIncreasing,
	      {G[ua,lb,lc],lb,lc}];
(* Print["Just after Map[AssertSymmetircDecreasing"]; *)
(******************************************************************************)
(* Shorthands *)
(******************************************************************************)
 
(* CD fixed error sizmaz -> sigmaz *)
q[rho1_,zcyl_] := -ampBrill*(rho1)^2 * ((rho1/sigmarho)^2 + (zcyl/sigmaz)^2);

(* 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 [g[la,lb]  ], prefix <> "g"],
   SetGroupName [CreateGroupFromTensor [gf[ua,ub]  ], prefix <> "gf"]
  };

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;
	
formgkCalc =  
{	
  Name -> "formgk",
  Schedule -> {"AT CCTK_INITIAL before CT_MultiLevel"},
  Where -> Interior, 
  Shorthands -> {rho1, e2q, rho2, detg, oneoverdetg},
  Equations ->
{
rho1 -> Sqrt[x*x + y*y],
e2q -> Exp[2*q[rho1,z]],
rho2 -> x*x + y*y,
g11 -> e2q + (one - e2q)*y*y/rho2,
g12 -> - (one - e2q)*x*y/rho2,
g13 -> 0,
g22 -> e2q + (one - e2q)*x*x/rho2,
g23 -> 0,
g33 -> e2q
}
};

 
detgExpr  = Det [MatrixOfComponents [g [la,lb]]];

gupperCalc =
{
  Name -> "gupper",
  Schedule -> {"AT CCTK_INITIAL before CT_MultiLevel AFTER formgkCalc"},
  Where -> Interior,
  Shorthands -> {},
  Equations ->
  {
    gf11 -> 1/g11 + g12 g12/(g11 g11 g22 -g11 g12 g12),
    gf12 -> -g12/(g11 g22 -g12 g12),
    gf13 -> 0,
    gf22 -> 1/(g22 - g12 g12 / g11),
    gf23 -> 0,
    gf33 -> 1/g33
  }
};

Print["Just before def of CT_BrillAnalyticCalc"];

BrillAnalyticCalc  =
{
  Name -> ThornName ,
  Schedule -> {"AT CCTK_INITIAL before CT_MultiLevel AFTER gupperCalc"},
(*  ConditionalOnKeyword -> {"free_data", "BrillWave_K"}, *)
  Shorthands -> {coeffpartialwrtxx, coeffpartialwrtxy, coeffpartialwrtyy, coeffpartialwrtzz, coeffpartialwrtx, coeffpartialwrty, coeffpartialwrtz, sqrtdetg, oneoversqrtdetg, detg, FF[ua], trRicci, Ricci[la,lb], gu[ua,ub], G[ua,lb,lc]},
  Where -> Interior,
  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 *)
    sqrtdetg -> Sqrt[detg],
    oneoversqrtdetg -> 1/sqrtdetg,
    FF[ua] -> oneoversqrtdetg ( PD[sqrtdetg gf[ua,ub], lb] ),
    coeffpartialwrtxx -> gf11,
    coeffpartialwrtxy -> gf12,
    coeffpartialwrtyy -> gf22,
    coeffpartialwrtzz -> gf33,
    coeffpartialwrtx  -> FF1,
    coeffpartialwrty  -> FF2,
    coeffpartialwrtz  -> FF3,
    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 *)
    Default -> 1/1000000
  }

};

(*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 =
{
  formgkCalc,
  gupperCalc,
  BrillAnalyticCalc
};

CreateKrancThornTT 
[
  groups, ".", ThornName,
  Calculations -> calculations,
  DeclaredGroups -> declaredGroupNames,
  PartialDerivatives -> derivatives,
  UseLoopControl -> True,
  UseVectors -> False,
  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];