# This is an until-end-of-line comment

THORN BSSN



DERIVATIVE
  PDstandardNth[i_]    -> StandardCenteredDifferenceOperator[1,derivOrder/2,i]

DERIVATIVE
  PDstandardNth[i_,i_] -> StandardCenteredDifferenceOperator[2,derivOrder/2,i]
DERIVATIVE
  PDstandardNth[i_,j_] -> StandardCenteredDifferenceOperator[1,derivOrder/2,i] *
                          StandardCenteredDifferenceOperator[1,derivOrder/2,j]
DERIVATIVE
  PDdissipationNth[i_] ->
    spacing[i]^(derivOrder+1) / 2^(derivOrder+2) *
    StandardCenteredDifferenceOperator[derivOrder+2,derivOrder/2+1,i]

JACOBIAN {PD, FD, J, dJ}



TENSOR normal, tangentA, tangentB, dir
TENSOR xx, rr, th, ph
TENSOR J, dJ
TENSOR admg, admK, admalpha, admdtalpha, admbeta, admdtbeta, H, M
TENSOR g, detg, gu, G, R, trR, Km, trK, cdphi, cdphi2
TENSOR phi, gt, At, Xt, Xtn, alpha, A, beta, B, Atm, Atu, trA, Ats, trAts
TENSOR dottrK, dotXt
TENSOR cXt, cS, cA
TENSOR e4phi, em4phi, ddetg, detgt, gtu, ddetgt, dgtu, ddgtu, Gtl, Gtlu, Gt
TENSOR Rt, Rphi, gK
TENSOR T00, T0, T, rho, S
TENSOR x, y, z, r
TENSOR epsdiss



SYMMETRIC admg[la,lb], admK[la,lb]
SYMMETRIC g[la,lb], K[la,lb], R[la,lb], cdphi2[la,lb]
SYMMETRIC gt[la,lb], At[la,lb], Ats[la,lb], Rt[la,lb], Rphi[la,lb], T[la,lb]

SYMMETRIC {dJ[ua,lb,lc], lb, lc}
SYMMETRIC {G[ua,lb,lc], lb, lc}
SYMMETRIC {Gtl[la,lb,lc], lb, lc}
SYMMETRIC {Gt[ua,lb,lc], lb, lc}
SYMMETRIC {gK[la,lb,lc], la, lb}

SYMMETRIC gu[ua,ub], gtu[ua,ub], Atu[ua,ub]
SYMMETRIC {dgtu[ua,ub,lc], ua, ub}
SYMMETRIC {ddgtu[ua,ub,lc,ld], ua, ub}
SYMMETRIC {ddgtu[ua,ub,lc,ld], lc, ld}



CONNECTION {CD, PD, G}
CONNECTION {CDt, PD, Gt}



GROUP {phi      , "log_confac"}
GROUP {gt[la,lb], "metric"    }
GROUP {Xt[ua   ], "Gamma"     }

EXTRA_GROUP {"Grid::coordinates", {x, y, z, r}}
EXTRA_GROUP {"ADMBase::metric",  {gxx, gxy, gxz, gyy, gyz, gzz}}
EXTRA_GROUP {"ADMBase::curv",    {kxx, kxy, kxz, kyy, kyz, kzz}}



DEFINE pi = N[Pi,40]

DEFINE
  detgExpr  = Det [MatrixOfComponents [g [la,lb]]]
DEFINE
  ddetgExpr[la_] =
    Sum [D[Det[MatrixOfComponents[g[la, lb]]], X] PD[X, la],
         {X, Union[Flatten[MatrixOfComponents[g[la, lb]]]]}]

DEFINE
  detgtExpr = Det [MatrixOfComponents [gt[la,lb]]]
DEFINE
  ddetgtExpr[la_] =
    Sum [D[Det[MatrixOfComponents[gt[la, lb]]], X] PD[X, la],
         {X, Union[Flatten[MatrixOfComponents[gt[la, lb]]]]}]



CALCULATION "Minkowski"
  Schedule: {"IN ADMBase_InitialData"}
  ConditionalOnKeyword: {"my_initial_data", "Minkowski"}
  BEGIN EQUATIONS
    phi       -> IfThen[conformalMethod, 1, 0]
    phi       -> conformalMethod ? 1 : 0
    gt[la,lb] -> KD[la,lb]
    trK       -> 0
    At[la,lb] -> 0
    Xt[ua]    -> 0
    alpha     -> 1
    A         -> 0
    beta[ua]  -> 0
    B[ua]     -> 0
  END EQUATIONS
END CALCULATION



CALCULATION "convertFromADMBase"
  Schedule: {"AT initial AFTER ADMBase_PostInitial"}
  ConditionalOnKeyword: {"my_initial_data", "ADMBase"}
  SHORTHAND g[la,lb], detg, gu[ua,ub], em4phi
  BEGIN EQUATIONS
    g[la,lb]  -> admg[la,lb]
    detg      -> detgExpr
    gu[ua,ub] -> 1/detg detgExpr MatrixInverse [g[ua,ub]]
    
    phi       -> IfThen [conformalMethod, detg^(-1/6), Log[detg]/12]
    em4phi    -> IfThen [conformalMethod, phi^2, Exp[-4 phi]]
    gt[la,lb] -> em4phi g[la,lb]
    
    trK       -> gu[ua,ub] admK[la,lb]
    At[la,lb] -> em4phi (admK[la,lb] - (1/3) g[la,lb] trK)
    
    alpha     -> admalpha
    
    beta[ua]  -> admbeta[ua]
  END EQUATIONS
END CALCULATION



CALCULATION convertFromADMBaseGammaCalc
  Name: BSSN <> "_convertFromADMBaseGamma"
  Schedule: {"AT initial AFTER " <> BSSN <> "_convertFromADMBase"}
  ConditionalOnKeyword: {"my_initial_data", "ADMBase"}
  # Do not synchronise right after this routine; instead, synchronise
  # after extrapolating
  Where: Interior
  # Synchronise after this routine, so that the refinement boundaries
  # are set correctly before extrapolating.  (We will need to
  # synchronise again after extrapolating because extrapolation does
  # not fill ghost zones, but this is irrelevant here.)
  SHORTHAND dir[ua]
  SHORTHAND detgt, gtu[ua,ub], Gt[ua,lb,lc], theta
  BEGIN EQUATIONS
    dir[ua] -> Sign[beta[ua]]
    
    detgt        -> 1 (* detgtExpr *)
    gtu[ua,ub]   -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]]
    Gt[ua,lb,lc] -> 1/2 gtu[ua,ud]
                    (PD[gt[lb,ld],lc] + PD[gt[lc,ld],lb] - PD[gt[lb,lc],ld])
    Xt[ua] -> gtu[ub,uc] Gt[ua,lb,lc]
    
    # If LapseACoeff=0, then A is not evolved, in the sense that it
    # does not influence the time evolution of other variables.
    A -> IfThen [LapseACoeff != 0
                 1 / (- harmonicF alpha^harmonicN)
                 (+ admdtalpha
                  - LapseAdvectionCoeff beta[ua] PDua[alpha,la]
                  - LapseAdvectionCoeff Abs[beta[ua]] PDus[alpha,la])
                 0]
    
    theta -> thetaExpr
    
    # If ShiftBCoeff=0 or theta ShiftGammaCoeff=0, then B^i is not
    # evolved, in the sense that it does not influence the time
    # evolution of other variables.
    B[ua] -> IfThen [ShiftGammaCoeff ShiftBCoeff != 0
                     1 / (theta ShiftGammaCoeff)
                     (+ admdtbeta[ua]
                      - ShiftAdvectionCoeff beta[ub] PDua[beta[ua],lb]
                      - ShiftAdvectionCoeff Abs[beta[ub]] PDus[beta[ua],lb])
                     0]
  END EQUATIONS
END CALCULATION



INHERITED_IMPLEMENTATION ADMBase, TmunuBase



KEYWORD_PARAMETER "my_initial_data"
  # Visibility: "restricted"
  # Description: "ddd"
  AllowedValues: {"ADMBase", "Minkowski"}
  Default: "ADMBase"
END KEYWORD_PARAMETER



REAL_PARAMETER LapseACoeff
  Description: "Whether to evolve A in time"
  Default: 0
END REAL_PARAMETER



END THORN