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];