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] }; 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, G, trRicci, gu }]; Map[AssertSymmetricDecreasing, { k[la,lb], g[la,lb] }]; (* Print["Just after Map[AssertSymmetircDecreasing"]; *) (******************************************************************************) (* Shorthands *) (******************************************************************************) q[rho1_,zcyl_] := -ampBrill*rho1*((rho1/sigmarho)^2 + (zcyl/sizmaz)^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 [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 *) (******************************************************************************) g11 = gxx; g12 = gxy; g13 = gxz; g22 = gyy; g23 = gyz; g33 = gzz; 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, Equations -> { 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; } }; TrRicciCalc = { Name -> "TrRicci", Schedule -> {"AT CCTK_INITIAL before CT_MultiLevel AFTER formgkCalc"}, 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 *) } }; Print["Just before def of CT_BrillAnalyticCalc"]; CT_BrillAnalyticCalc = { Name -> ThornName , Schedule -> {"AT CCTK_INITIAL before CT_MultiLevel AFTER TrRicciCalc"}, (* ConditionalOnKeyword -> {"free_data", "BrillWave_K"}, *) Where -> Interior, 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 *) 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, TrRicciCalc, CT_BrillAnalyticCalc }; 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];