From 40ee6d1f89ecd9bd370a5310cfc7bb50a2b81910 Mon Sep 17 00:00:00 2001 From: Alex Johnson Date: Wed, 28 Jun 2017 17:32:16 +0530 Subject: [PATCH 1/3] OTWO-3061 Fix mathematica comment parsing The original code in PR#23 did not consider multiline comments. We compared and copied the code from OCaml parser which has a similar comment structure and it worked. The existing test files in test/expected_dir & test/src_dir needed fixing for some mismatch in whitespace characters and expected lines. --- src/parsers/mathematica.rl | 36 +- test/expected_dir/mathematica1.m | 118 ++- test/src_dir/mathematica1.m | 1480 ++++++++++++++++++++++++++++++++++++++ 3 files changed, 1564 insertions(+), 70 deletions(-) create mode 100644 test/src_dir/mathematica1.m diff --git a/src/parsers/mathematica.rl b/src/parsers/mathematica.rl index dcf142e..72b5957 100644 --- a/src/parsers/mathematica.rl +++ b/src/parsers/mathematica.rl @@ -44,14 +44,24 @@ enum { } } - mathematica_comment = - '(*' @comment ( - newline %{ entity = INTERNAL_NL; } %mathematica_ccallback - | - ws - | - (nonnewline - ws) @code - )* :>> '*)'; + action mathematica_comment_nc_res { nest_count = 0; } + action mathematica_comment_nc_inc { nest_count++; } + action mathematica_comment_nc_dec { nest_count--; } + + mathematica_nested_block_comment = + '(*' >mathematica_comment_nc_res @comment ( + newline %{ entity = INTERNAL_NL; } %mathematica_ccallback + | + ws + | + '(*' @mathematica_comment_nc_inc @comment + | + '*)' @mathematica_comment_nc_dec @comment + | + (nonnewline - ws) @comment + )* :>> ('*)' when { nest_count == 0 }) @comment; + + mathematica_comment = mathematica_nested_block_comment; mathematica_string = '"' @code ([^"]) '"'; @@ -70,7 +80,13 @@ enum { userdata); } - mathematica_comment_entity = '(*' any* :>> '*)'; + mathematica_comment_entity = '(*' >mathematica_comment_nc_res ( + '(*' @mathematica_comment_nc_inc + | + '*)' @mathematica_comment_nc_dec + | + any + )* :>> ('*)' when { nest_count == 0 }); mathematica_entity := |* space+ ${ entity = MATHEMATICA_SPACE; } => mathematica_ecallback; @@ -100,6 +116,8 @@ void parse_mathematica(char *buffer, int length, int count, ) { init + int nest_count = 0; + %% write init; cs = (count) ? mathematica_en_mathematica_line : mathematica_en_mathematica_entity; %% write exec; diff --git a/test/expected_dir/mathematica1.m b/test/expected_dir/mathematica1.m index b48ed57..43fcb2a 100644 --- a/test/expected_dir/mathematica1.m +++ b/test/expected_dir/mathematica1.m @@ -42,7 +42,7 @@ mathematica comment (*********************************************************** mathematica code (-1)^(fdOrder/2) * mathematica code spacing[i]^(fdOrder+1) / 2^(fdOrder+2) * mathematica code StandardCenteredDifferenceOperator[fdOrder+2,fdOrder/2+1,i], -mathematica blank +mathematica blank mathematica comment (* PD: These come from my mathematica notebook mathematica comment "Upwind-Kranc-Convert.nb" that converts upwinding finite mathematica comment differencing operators generated by @@ -294,16 +294,16 @@ mathematica comment (*********************************************************** mathematica code g[la,lb] -> admg[la,lb], mathematica code detg -> detgExpr, mathematica code gu[ua,ub] -> 1/detg detgExpr MatrixInverse [g[ua,ub]], -mathematica code +mathematica blank mathematica code phi -> IfThen[conformalMethod==CMW, detg^(-1/6), Log[detg]/12], mathematica code em4phi -> IfThen[conformalMethod==CMW, phi^2, Exp[-4 phi]], mathematica code gt[la,lb] -> em4phi g[la,lb], -mathematica code +mathematica blank mathematica code trK -> gu[ua,ub] admK[la,lb], mathematica code At[la,lb] -> em4phi (admK[la,lb] - (1/3) g[la,lb] trK), -mathematica code +mathematica blank mathematica code alpha -> admalpha, -mathematica code +mathematica blank mathematica code beta[ua] -> admbeta[ua], mathematica blank mathematica code IfCCZ4[Theta -> 0] @@ -330,13 +330,13 @@ mathematica comment (* Synchronise after this routine, so that the refinement mathematica code Equations -> mathematica code { mathematica code dir[ua] -> Sign[beta[ua]], -mathematica code +mathematica blank mathematica code detgt -> 1 (* detgtExpr *), mathematica code gtu[ua,ub] -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]], mathematica code Gt[ua,lb,lc] -> 1/2 gtu[ua,ud] mathematica code (PD[gt[lb,ld],lc] + PD[gt[lc,ld],lb] - PD[gt[lb,lc],ld]), mathematica code Xt[ua] -> gtu[ub,uc] Gt[ua,lb,lc], -mathematica code +mathematica blank mathematica comment (* mathematica comment A -> - admdtalpha / (harmonicF alpha^harmonicN) (LapseAdvectionCoeff - 1), mathematica comment *) @@ -347,9 +347,9 @@ mathematica comment (* If LapseACoeff=0, then A is not evolved, in the sense mathematica code (+ admdtalpha mathematica code - LapseAdvectionCoeff Upwind[beta[ua], alpha, la]), mathematica code 0], -mathematica code +mathematica blank mathematica code theta -> thetaExpr, -mathematica code +mathematica blank mathematica comment (* If ShiftBCoeff=0 or theta ShiftGammaCoeff=0, then B^i is not mathematica comment evolved, in the sense that it does not influence the time mathematica comment evolution of other variables. *) @@ -413,14 +413,14 @@ mathematica comment (*********************************************************** mathematica code Equations -> mathematica code { mathematica code dir[ua] -> Sign[beta[ua]], -mathematica code +mathematica blank mathematica code detgt -> 1 (* detgtExpr *), mathematica comment (* This leads to simpler code... *) mathematica code gtu[ua,ub] -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]], -mathematica code +mathematica blank mathematica code eta -> etaExpr, mathematica code theta -> thetaExpr, -mathematica code +mathematica blank mathematica comment (* see RHS *) mathematica comment (* mathematica comment admdtalpha -> - harmonicF alpha^harmonicN @@ -438,7 +438,7 @@ mathematica code (- 2 alpha PD[phi,lj] mathematica code + 2 phi PD[alpha,lj] mathematica code + gtu[uk,ul] phi alpha mathematica code (PD[gt[lk,ll],lj] - 2 PD[gt[lj,lk],ll])), -mathematica code (* else *) +mathematica comment (* else *) mathematica code + theta ShiftGammaCoeff mathematica code (+ ShiftBCoeff B[ua] mathematica code + (1 - ShiftBCoeff) @@ -459,10 +459,10 @@ mathematica code (Xt[ua] - eta BetaDriver beta[ mathematica code detgt -> 1 (* detgtExpr *), mathematica comment (* This leads to simpler code... *) mathematica code gtu[ua,ub] -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]], -mathematica code +mathematica blank mathematica code eta -> etaExpr, mathematica code theta -> thetaExpr, -mathematica code +mathematica blank mathematica comment (* see RHS, but omit derivatives near the boundary *) mathematica comment (* mathematica comment admdtalpha -> - harmonicF alpha^harmonicN @@ -474,7 +474,7 @@ mathematica code (+ LapseACoeff A mathematica code (trK - IfCCZ4[2 Theta, 0]))), mathematica code admdtbeta[ua] -> IfThen[harmonicShift, mathematica code 0, -mathematica code (* else *) +mathematica comment (* else *) mathematica code + theta ShiftGammaCoeff mathematica code (+ ShiftBCoeff B[ua] mathematica code + (1 - ShiftBCoeff) @@ -494,10 +494,10 @@ mathematica code (Xt[ua] - eta BetaDriver beta[ mathematica code detgt -> 1 (* detgtExpr *), mathematica comment (* This leads to simpler code... *) mathematica code gtu[ua,ub] -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]], -mathematica code +mathematica blank mathematica code eta -> etaExpr, mathematica code theta -> thetaExpr, -mathematica code +mathematica blank mathematica comment (* see RHS, but omit derivatives everywhere (which is wrong, but mathematica comment faster, since it does not require synchronisation or boundary mathematica comment conditions) *) @@ -511,7 +511,7 @@ mathematica code (+ LapseACoeff A mathematica code (trK - IfCCZ4[2 Theta, 0]))), mathematica code admdtbeta[ua] -> IfThen[harmonicShift, mathematica code 0, -mathematica code (* else *) +mathematica comment (* else *) mathematica code + theta ShiftGammaCoeff mathematica code (+ ShiftBCoeff B[ua] mathematica code + (1 - ShiftBCoeff) @@ -547,16 +547,16 @@ mathematica comment (* Synchronise the RHS grid functions after this routine, mathematica code Equations -> mathematica code { mathematica code dir[ua] -> Sign[beta[ua]], -mathematica code +mathematica blank mathematica code detgt -> 1 (* detgtExpr *), -mathematica code +mathematica blank mathematica comment (* This leads to simpler code... *) mathematica code gtu[ua,ub] -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]], mathematica code Gtl[la,lb,lc] -> 1/2 mathematica code (PD[gt[lb,la],lc] + PD[gt[lc,la],lb] - PD[gt[lb,lc],la]), mathematica code Gtlu[la,lb,uc] -> gtu[uc,ud] Gtl[la,lb,ld], mathematica code Gt[ua,lb,lc] -> gtu[ua,ud] Gtl[ld,lb,lc], -mathematica code +mathematica blank mathematica comment (* The conformal connection functions calculated from the conformal metric, mathematica comment used instead of Xt where no derivatives of Xt are taken *) mathematica code Xtn[ui] -> gtu[uj,uk] Gt[ui,lj,lk], @@ -594,36 +594,36 @@ mathematica comment (* PRD 62, 044034 (2000), eqn. (15) *) mathematica code - 2 gt[li,lj] gtu[ul,un] cdphi2[ll,ln] mathematica code + 4 cdphi[li] cdphi[lj] mathematica code - 4 gt[li,lj] gtu[ul,un] cdphi[ln] cdphi[ll], -mathematica code +mathematica blank mathematica code Atm[ua,lb] -> gtu[ua,uc] At[lc,lb], mathematica code Atu[ua,ub] -> Atm[ua,lc] gtu[ub,uc], -mathematica code +mathematica blank mathematica code R[la,lb] -> Rt[la,lb] + Rphi[la,lb], mathematica code IfCCZ4[ mathematica code R[la,lb] -> R[la,lb] + (2/phi) (+ g[la,lc] Z[uc] PD[phi,lb] mathematica code + g[lb,lc] Z[uc] PD[phi,la] - g[la,lb] Z[uc] PD[phi,lc]) mathematica code + e4phi Z[uc] PD[gt[la,lb],lc] mathematica code ], -mathematica code +mathematica blank mathematica comment (* Matter terms *) -mathematica blank +mathematica blank mathematica comment (* rho = n^a n^b T_ab *) mathematica code rho -> addMatter mathematica code (1/alpha^2 (T00 - 2 beta[ui] T0[li] + beta[ui] beta[uj] T[li,lj])), -mathematica code +mathematica blank mathematica comment (* S_i = -p^a_i n^b T_ab, where p^a_i = delta^a_i + n^a n_i *) mathematica code S[li] -> addMatter (-1/alpha (T0[li] - beta[uj] T[li,lj])), -mathematica code +mathematica blank mathematica comment (* trS = gamma^ij T_ij *) mathematica code trS -> addMatter (em4phi gtu[ui,uj] T[li,lj]), -mathematica code +mathematica blank mathematica comment (* RHS terms *) -mathematica code +mathematica blank mathematica comment (* PRD 62, 044034 (2000), eqn. (10) *) mathematica comment (* PRD 67 084023 (2003), eqn. (16) and (23) *) mathematica code dot[phi] -> IfThen[conformalMethod==CMW, 1/3 phi, -1/6] mathematica code (alpha trK - PD[beta[ua],la]), -mathematica code +mathematica blank mathematica comment (* PRD 62, 044034 (2000), eqn. (9) *) mathematica comment (* gr-qc:1106.2254 (2011), eqn. (14) *) mathematica comment (* removing trA from Aij ensures that detg = 1 *) @@ -700,17 +700,14 @@ mathematica comment (* Adding Z terms in the Ricci and Theta terms *) mathematica comment (* Equation (4.23) in Baumgarte & Shapiro (Phys. Rept. 376 (2003) 41-131) *) mathematica code + addMatter (- em4phi alpha 8 Pi mathematica code (T[la,lb] - (1/3) g[la,lb] trS)), -mathematica code +mathematica blank mathematica comment (* dot[alpha] -> - harmonicF alpha^harmonicN trK, *) mathematica comment (* dot[alpha] -> - harmonicF alpha^harmonicN A + Lie[alpha, beta], *) mathematica comment (* mathematica comment dot[alpha] -> - harmonicF alpha^harmonicN ( mathematica comment (1 - LapseAdvectionCoeff) A + LapseAdvectionCoeff trK) mathematica comment + LapseAdvectionCoeff beta[ua] PDu[alpha,la], -mathematica blank -mathematica blank mathematica comment dot[A] -> (1 - LapseAdvectionCoeff) (dottrK - AlphaDriver A), -mathematica blank mathematica comment *) mathematica code dot[alpha] -> - harmonicF alpha^harmonicN mathematica code (+ LapseACoeff A @@ -721,10 +718,10 @@ mathematica code (+ trK - IfCCZ4[2 Theta, 0] mathematica code dot[A] -> + (LapseACoeff mathematica code (+ dottrK - IfCCZ4[2 dotTheta, 0] mathematica code - AlphaDriver A)), -mathematica code +mathematica blank mathematica code eta -> etaExpr, mathematica code theta -> thetaExpr, -mathematica code +mathematica blank mathematica comment (* dot[beta[ua]] -> eta Xt[ua], *) mathematica comment (* dot[beta[ua]] -> ShiftGammaCoeff alpha^ShiftAlphaPower B[ua], *) mathematica code dot[beta[ua]] -> IfThen[harmonicShift, @@ -939,20 +936,20 @@ mathematica comment (* Initialise the RHS variables in analysis in case they are mathematica code Equations -> mathematica code { mathematica code dir[ua] -> Sign[normal[ua]], -mathematica code +mathematica blank mathematica code detgt -> 1 (* detgtExpr *), mathematica code gtu[ua,ub] -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]], mathematica code em4phi -> IfThen[conformalMethod==CMW, phi^2, Exp[-4 phi]], mathematica code gu[ua,ub] -> em4phi gtu[ua,ub], -mathematica code +mathematica blank mathematica code nn[la] -> Euc[la,lb] normal[ub], mathematica code nu[ua] -> gu[ua,ub] nn[lb], mathematica code nlen2 -> nu[ua] nn[la], mathematica code nlen -> Sqrt[nlen2], mathematica code su[ua] -> nu[ua] / nlen, -mathematica code +mathematica blank mathematica code vg -> Sqrt[harmonicF], -mathematica code +mathematica blank mathematica code dot[phi] -> - vg su[uc] PDo[phi ,lc], mathematica code dot[gt[la,lb]] -> - su[uc] PDo[gt[la,lb],lc], mathematica code dot[trK] -> - vg su[uc] PDo[trK ,lc], @@ -986,11 +983,11 @@ mathematica comment which is not always the case (since we don't enforce mathematica comment the other hand, this may not be so important... *) mathematica code detgt -> 1 (* detgtExpr *), mathematica code gtu[ua,ub] -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]], -mathematica code +mathematica blank mathematica code trAt -> gtu[ua,ub] At[la,lb], -mathematica code +mathematica blank mathematica code At[la,lb] -> At[la,lb] - (1/3) gt[la,lb] trAt, -mathematica code +mathematica blank mathematica code alpha -> Max[alpha, MinimumLapse] mathematica code } mathematica code }; @@ -1041,14 +1038,14 @@ mathematica comment (*********************************************************** mathematica code { mathematica code detgt -> 1 (* detgtExpr *), mathematica code ddetgt[la] -> 0 (* ddetgtExpr[la] *), -mathematica code +mathematica blank mathematica comment (* This leads to simpler code... *) mathematica code gtu[ua,ub] -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]], mathematica code Gtl[la,lb,lc] -> 1/2 mathematica code (PD[gt[lb,la],lc] + PD[gt[lc,la],lb] - PD[gt[lb,lc],la]), mathematica code Gtlu[la,lb,uc] -> gtu[uc,ud] Gtl[la,lb,ld], mathematica code Gt[ua,lb,lc] -> gtu[ua,ud] Gtl[ld,lb,lc], -mathematica code +mathematica blank mathematica comment (* The conformal connection functions calculated from the conformal metric, mathematica comment used instead of Xt where no derivatives of Xt are taken *) mathematica code Xtn[ui] -> gtu[uj,uk] Gt[ui,lj,lk], @@ -1063,7 +1060,7 @@ mathematica comment (* The Z quantities *) mathematica code IfCCZ4[ mathematica code Z[ud] -> (1/2) gu[ua,ud] (- PD[gt[la,lb],lc] gtu[ub,uc] + gt[la,lc] Xt[uc]) mathematica code ], -mathematica code +mathematica blank mathematica comment (* PRD 62, 044034 (2000), eqn. (18) *) mathematica code Rt[li,lj] -> - (1/2) gtu[ul,um] PD[gt[li,lj],ll,lm] mathematica code + (1/2) gt[lk,li] PD[Xt[uk],lj] @@ -1106,14 +1103,14 @@ mathematica comment (* PRD 62, 044034 (2000), eqn. (15) *) mathematica code - 2 gt[li,lj] gtu[ul,un] cdphi2[ll,ln] mathematica code + 4 cdphi[li] cdphi[lj] mathematica code - 4 gt[li,lj] gtu[ul,un] cdphi[ln] cdphi[ll], -mathematica code +mathematica blank mathematica comment (* ddetg[la] -> PD[e4phi detg,la], *) mathematica code ddetg[la] -> e4phi ddetgt[la] + 4 detgt e4phi PD[phi,la], mathematica comment (* TODO: check this equation, maybe simplify it by omitting ddetg *) mathematica code G[ua,lb,lc] -> Gt[ua,lb,lc] mathematica code + 1/(2 detg) (+ KD[ua,lb] ddetg[lc] + KD[ua,lc] ddetg[lb] mathematica code - (1/3) g[lb,lc] gu[ua,ud] ddetg[ld]), -mathematica code +mathematica blank mathematica code R[la,lb] -> + Rt[la,lb] + Rphi[la,lb], mathematica blank mathematica code IfCCZ4[ @@ -1123,42 +1120,41 @@ mathematica comment (* TODO: check this equation, maybe simplify it by omitt mathematica code ], mathematica blank mathematica code trR -> gu[ua,ub] R[la,lb], -mathematica code +mathematica blank mathematica comment (* K[la,lb] -> e4phi At[la,lb] + (1/3) g[la,lb] trK, *) mathematica comment (* Km[ua,lb] -> gu[ua,uc] K[lc,lb], *) mathematica code Atm[ua,lb] -> gtu[ua,uc] At[lc,lb], -mathematica code +mathematica blank mathematica comment (* Matter terms *) -mathematica code +mathematica blank mathematica comment (* rho = n^a n^b T_ab *) mathematica code rho -> 1/alpha^2 (T00 - 2 beta[ui] T0[li] + beta[ui] beta[uj] T[li,lj]), -mathematica code +mathematica blank mathematica comment (* S_i = -p^a_i n^b T_ab, where p^a_i = delta^a_i + n^a n_i *) mathematica code S[li] -> -1/alpha (T0[li] - beta[uj] T[li,lj]), -mathematica code +mathematica blank mathematica comment (* Constraints *) -mathematica code +mathematica blank mathematica comment (* H -> trR - Km[ua,lb] Km[ub,la] + trK^2, *) mathematica comment (* PRD 67, 084023 (2003), eqn. (19) *) mathematica code H -> trR - Atm[ua,lb] Atm[ub,la] + (2/3) trK^2 - addMatter 16 Pi rho, -mathematica code +mathematica blank mathematica comment (* gK[la,lb,lc] -> CD[K[la,lb],lc], *) mathematica comment (* gK[la,lb,lc] -> + 4 e4phi PD[phi,lc] At[la,lb] + e4phi CD[At[la,lb],lc] mathematica comment + (1/3) g[la,lb] PD[trK,lc], -mathematica blank mathematica comment M[la] -> gu[ub,uc] (gK[lc,la,lb] - gK[lc,lb,la]), *) mathematica blank mathematica code M[li] -> + gtu[uj,uk] (CDt[At[li,lj],lk] + 6 At[li,lj] cdphi[lk]) mathematica code - (2/3) PD[trK,li] mathematica code - addMatter 8 Pi S[li], mathematica comment (* TODO: use PRD 67, 084023 (2003), eqn. (20) *) -mathematica code +mathematica blank mathematica comment (* det gamma-tilde *) mathematica code cS -> Log[detgt], -mathematica code +mathematica blank mathematica comment (* Gamma constraint *) mathematica code cXt[ua] -> gtu[ub,uc] Gt[ua,lb,lc] - Xt[ua], -mathematica code +mathematica blank mathematica comment (* trace A-tilde *) mathematica code cA -> gtu[ua,ub] At[la,lb] mathematica code } @@ -1481,4 +1477,4 @@ mathematica comment (* These are the arguments to createComment: mathematica code createCode[4, False, True, True , True , 3, 1, "BSSN"]; mathematica blank mathematica code createCode[4, False, True, True , False, 3, 1, "CCZ4"]; -mathematica blank \ No newline at end of file +mathematica blank diff --git a/test/src_dir/mathematica1.m b/test/src_dir/mathematica1.m new file mode 100644 index 0000000..4bec755 --- /dev/null +++ b/test/src_dir/mathematica1.m @@ -0,0 +1,1480 @@ + +SetEnhancedTimes[False]; +SetSourceLanguage["C"]; + +(******************************************************************************) +(* Options *) +(******************************************************************************) + +createCode[derivOrder_, useJacobian_, splitUpwindDerivs_, useVectors_, useOpenCL_, evolutionTimelevels_, addMatter_, formulation_] := +Module[{prefix, suffix, thorn}, + +prefix = "ML_"; +suffix = + "" + <> If [useJacobian, "_MP", ""] + <> If [derivOrder!=4, "_O" <> ToString[derivOrder], ""] + <> If [splitUpwindDerivs, "", "_UPW"] + <> If [useVectors, "", "_NV"] + <> If [useOpenCL, "_CL", ""] + (* <> If [evolutionTimelevels!=3, "_TL" <> ToString[evolutionTimelevels], ""] *) + (* <> If [addMatter==1, "_M", ""] *) + ; + +thorn = prefix <> formulation <> suffix; + +SetAttributes[IfCCZ4, HoldAll]; +IfCCZ4[expr_, else_:Sequence[]] := If[formulation === "CCZ4", expr, Unevaluated[else]]; + +(******************************************************************************) +(* Derivatives *) +(******************************************************************************) + +KD = KroneckerDelta; + +derivatives = +{ + PDstandardNth[i_] -> StandardCenteredDifferenceOperator[1,fdOrder/2,i], + PDstandardNth[i_,i_] -> StandardCenteredDifferenceOperator[2,fdOrder/2,i], + PDstandardNth[i_,j_] -> StandardCenteredDifferenceOperator[1,fdOrder/2,i] * + StandardCenteredDifferenceOperator[1,fdOrder/2,j], + PDdissipationNth[i_] -> + (-1)^(fdOrder/2) * + spacing[i]^(fdOrder+1) / 2^(fdOrder+2) * + StandardCenteredDifferenceOperator[fdOrder+2,fdOrder/2+1,i], + +(* PD: These come from my mathematica notebook + "Upwind-Kranc-Convert.nb" that converts upwinding finite + differencing operators generated by + StandardUpwindDifferenceOperator into this form *) + + Sequence@@Flatten[Table[ + {PDupwindNth[i] -> Switch[fdOrder, + 2, (dir[i]*(-3 + 4*shift[i]^dir[i] - shift[i]^(2*dir[i])))/(2*spacing[i]), + 4, (dir[i]*(-10 - 3/shift[i]^dir[i] + 18*shift[i]^dir[i] - + 6*shift[i]^(2*dir[i]) + shift[i]^(3*dir[i])))/(12*spacing[i]), + 6, (dir[i]*(-35 + 2/shift[i]^(2*dir[i]) - 24/shift[i]^dir[i] + 80*shift[i]^dir[i] - + 30*shift[i]^(2*dir[i]) + 8*shift[i]^(3*dir[i]) - shift[i]^(4*dir[i])))/(60*spacing[i]), + 8, (dir[i]*(-378 - 5/shift[i]^(3*dir[i]) + 60/shift[i]^(2*dir[i]) - 420/shift[i]^dir[i] + + 1050*shift[i]^dir[i] - 420*shift[i]^(2*dir[i]) + 140*shift[i]^(3*dir[i]) - 30*shift[i]^(4*dir[i]) + + 3*shift[i]^(5*dir[i])))/(840*spacing[i])], + + PDupwindNthAnti[i] -> Switch[fdOrder, + 2, (+1 shift[i]^(-2) -4 shift[i]^(-1) +0 shift[i]^( 0) +4 shift[i]^(+1) -1 shift[i]^(+2)) / (4 spacing[i]), + 4, (-1 shift[i]^(-3) +6 shift[i]^(-2) -21 shift[i]^(-1 )+0 shift[i]^( 0) +21 shift[i]^(+1) + -6 shift[i]^(+2) +1 shift[i]^(+3)) / (24 spacing[i]), + 6, (+1 shift[i]^(-4) -8 shift[i]^(-3) +32 shift[i]^(-2) -104 shift[i]^(-1) +0 shift[i]^( 0) + +104 shift[i]^(+1) -32 shift[i]^(+2) +8 shift[i]^(+3) -1 shift[i]^(+4)) / (120 spacing[i]), + 8, (-3 shift[i]^(-5) +30 shift[i]^(-4) -145 shift[i]^(-3) +480 shift[i]^(-2) -1470 shift[i]^(-1) + +0 shift[i]^( 0) +1470 shift[i]^(+1) -480 shift[i]^(+2) +145 shift[i]^(+3) -30 shift[i]^(+4) + +3 shift[i]^(+5)) / (1680 spacing[i])], + + PDupwindNthSymm[i] -> Switch[fdOrder, + 2, (-1 shift[i]^(-2) +4 shift[i]^(-1) -6 shift[i]^( 0) +4 shift[i]^(+1) -1 shift[i]^(+2)) / (4 spacing[i]), + 4, (+1 shift[i]^(-3) -6 shift[i]^(-2) +15 shift[i]^(-1) -20 shift[i]^( 0) +15 shift[i]^(+1) + -6 shift[i]^(+2) +1 shift[i]^(+3)) / (24 spacing[i]), + 6, (-1 shift[i]^(-4) +8 shift[i]^(-3) - 28 shift[i]^(-2)+56 shift[i]^(-1)-70 shift[i]^( 0) + +56 shift[i]^(+1) -28 shift[i]^(+2) +8 shift[i]^(+3) -1 shift[i]^(+4)) / (120 spacing[i]), + 8, (+3 shift[i]^(-5) -30 shift[i]^(-4) +135 shift[i]^(-3) -360 shift[i]^(-2) +630 shift[i]^(-1) + -756 shift[i]^( 0) +630 shift[i]^(+1) -360 shift[i]^(+2) +135 shift[i]^(+3) -30 shift[i]^(+4) + +3 shift[i]^(+5)) / (1680 spacing[i])], + + (* TODO: make these higher order stencils *) + PDonesided[i] -> dir[i] (-1 + shift[i]^dir[i]) / spacing[i]} /. i->j, {j,1,3}],1] +}; + +PD = PDstandardNth; +PDu = PDupwindNth; +PDua = PDupwindNthAnti; +PDus = PDupwindNthSymm; +(* PDo = PDonesided; *) +PDdiss = PDdissipationNth; + +If [splitUpwindDerivs, + Upwind[dir_, var_, idx_] := dir PDua[var,idx] + Abs[dir] PDus[var,idx], + Upwind[dir_, var_, idx_] := dir PDu[var,idx]]; + + + +(******************************************************************************) +(* Tensors *) +(******************************************************************************) + +(* Register the tensor quantities with the TensorTools package *) +Map [DefineTensor, + {normal, tangentA, tangentB, dir, + nn, nu, nlen, nlen2, su, vg, + xx, rr, th, ph, + admg, admK, admalpha, admdtalpha, admbeta, admdtbeta, H, M, + g, detg, gu, G, R, trR, Km, trK, cdphi, cdphi2, + phi, gt, At, Xt, Xtn, Theta, Z, + alpha, A, beta, B, Atm, Atu, trA, Ats, trAts, + dottrK, dotXt, + cXt, cS, cA, + e4phi, em4phi, ddetg, detgt, gtu, ddetgt, dgtu, ddgtu, Gtl, Gtlu, Gt, + Rt, Rphi, gK, + T00, T0, T, rho, S, + x, y, z, r, + epsdiss}]; + +(* NOTE: It seems as if Lie[.,.] did not take these tensor weights + into account. Presumably, CD[.,.] and CDt[.,.] don't do this either. *) +SetTensorAttribute[phi, TensorWeight, +1/6]; +SetTensorAttribute[gt, TensorWeight, -2/3]; +SetTensorAttribute[Xt, TensorWeight, +2/3]; +SetTensorAttribute[At, TensorWeight, -2/3]; +SetTensorAttribute[cXt, TensorWeight, +2/3]; +SetTensorAttribute[cS, TensorWeight, +2 ]; + +Map [AssertSymmetricIncreasing, + {admg[la,lb], admK[la,lb], g[la,lb], K[la,lb], R[la,lb], cdphi2[la,lb], + gt[la,lb], At[la,lb], Ats[la,lb], Rt[la,lb], Rphi[la,lb], T[la,lb]}]; +AssertSymmetricIncreasing [G[ua,lb,lc], lb, lc]; +AssertSymmetricIncreasing [Gtl[la,lb,lc], lb, lc]; +AssertSymmetricIncreasing [Gt[ua,lb,lc], lb, lc]; +AssertSymmetricIncreasing [gK[la,lb,lc], la, lb]; +Map [AssertSymmetricIncreasing, + {gu[ua,ub], gtu[ua,ub], Atu[ua,ub]}]; +AssertSymmetricIncreasing [dgtu[ua,ub,lc], ua, ub]; +AssertSymmetricIncreasing [ddgtu[ua,ub,lc,ld], ua, ub]; +AssertSymmetricIncreasing [ddgtu[ua,ub,lc,ld], lc, ld]; + +DefineConnection [CD, PD, G]; +DefineConnection [CDt, PD, Gt]; + +(* Use the CartGrid3D variable names *) +x1=x; x2=y; x3=z; + +(* Use the ADMBase variable names *) +admg11=gxx; admg12=gxy; admg22=gyy; admg13=gxz; admg23=gyz; admg33=gzz; +admK11=kxx; admK12=kxy; admK22=kyy; admK13=kxz; admK23=kyz; admK33=kzz; +admalpha=alp; +admdtalpha=dtalp; +admbeta1=betax; admbeta2=betay; admbeta3=betaz; +admdtbeta1=dtbetax; admdtbeta2=dtbetay; admdtbeta3=dtbetaz; + +(* Use the TmunuBase variable names *) +T00=eTtt; +T01=eTtx; T02=eTty; T03=eTtz; +T11=eTxx; T12=eTxy; T22=eTyy; T13=eTxz; T23=eTyz; T33=eTzz; + + + +(******************************************************************************) +(* Expressions *) +(******************************************************************************) + +(* enum constants for conformalMethod; these must be consistent + with the definition of the Cactus parameter conformalMethod *) +CMphi = 0; +CMW = 1; + +detgExpr = Det [MatrixOfComponents [g [la,lb]]]; +ddetgExpr[la_] = + Sum [D[Det[MatrixOfComponents[g[la, lb]]], X] PD[X, la], + {X, Union[Flatten[MatrixOfComponents[g[la, lb]]]]}]; + +detgtExpr = Det [MatrixOfComponents [gt[la,lb]]]; +ddetgtExpr[la_] = + Sum [D[Det[MatrixOfComponents[gt[la, lb]]], X] PD[X, la], + {X, Union[Flatten[MatrixOfComponents[gt[la, lb]]]]}]; + +etaExpr = SpatialBetaDriverRadius / Max [r, SpatialBetaDriverRadius]; +thetaExpr = Min [Exp [1 - r / SpatialShiftGammaCoeffRadius], 1]; + + + +(******************************************************************************) +(* Groups *) +(******************************************************************************) + +evolvedGroups = + {SetGroupName [CreateGroupFromTensor [phi ], prefix <> "log_confac"], + SetGroupName [CreateGroupFromTensor [gt[la,lb]], prefix <> "metric" ], + SetGroupName [CreateGroupFromTensor [Xt[ua] ], prefix <> "Gamma" ], + SetGroupName [CreateGroupFromTensor [trK ], prefix <> "trace_curv"], + SetGroupName [CreateGroupFromTensor [At[la,lb]], prefix <> "curv" ], + SetGroupName [CreateGroupFromTensor [alpha ], prefix <> "lapse" ], + SetGroupName [CreateGroupFromTensor [A ], prefix <> "dtlapse" ], + SetGroupName [CreateGroupFromTensor [beta[ua] ], prefix <> "shift" ], + SetGroupName [CreateGroupFromTensor [B[ua] ], prefix <> "dtshift" ], + IfCCZ4[SetGroupName[CreateGroupFromTensor[Theta], prefix <> "Theta"]]}; +evaluatedGroups = + {SetGroupName [CreateGroupFromTensor [H ], prefix <> "Ham"], + SetGroupName [CreateGroupFromTensor [M[la] ], prefix <> "mom"], + SetGroupName [CreateGroupFromTensor [cS ], prefix <> "cons_detg"], + SetGroupName [CreateGroupFromTensor [cXt[ua]], prefix <> "cons_Gamma"], + SetGroupName [CreateGroupFromTensor [cA ], prefix <> "cons_traceA"]}; + +declaredGroups = Join [evolvedGroups, evaluatedGroups]; +declaredGroupNames = Map [First, declaredGroups]; + + + +extraGroups = + {{"grid::coordinates", {x, y, z, r}}, + {"ADMBase::metric", {gxx, gxy, gxz, gyy, gyz, gzz}}, + {"ADMBase::curv", {kxx, kxy, kxz, kyy, kyz, kzz}}, + {"ADMBase::lapse", {alp}}, + {"ADMBase::dtlapse", {dtalp}}, + {"ADMBase::shift", {betax, betay, betaz}}, + {"ADMBase::dtshift", {dtbetax, dtbetay, dtbetaz}}, + {"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]; + + + +(******************************************************************************) +(* Initial data *) +(******************************************************************************) + +initialCalc = +{ + Name -> thorn <> "_Minkowski", + Schedule -> {"IN ADMBase_InitialData"}, + ConditionalOnKeyword -> {"my_initial_data", "Minkowski"}, + Equations -> + { + phi -> IfThen[conformalMethod==CMW, 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, + IfCCZ4[Theta -> 0] + } +}; + + + +(******************************************************************************) +(* Split a calculation *) +(******************************************************************************) + +PartialCalculation[calc_, suffix_, updates_, evolVars_] := +Module[ + {name, calc1, replaces, calc2, vars, patterns, eqs, calc3}, + (* Add suffix to name *) + name = lookup[calc, Name] <> suffix; + calc1 = mapReplace[calc, Name, name]; + (* Replace some entries in the calculation *) + (* replaces = Map[Function[rule, mapReplace[#, rule[[1]], rule[[2]]]&], updates]; *) + replaces = updates //. (lhs_ -> rhs_) -> (mapReplace[#, lhs, rhs]&); + calc2 = Apply[Composition, replaces][calc1]; + (* Remove unnecessary equations *) + vars = Join[evolVars, lookup[calc2, Shorthands]]; + patterns = Replace[vars, { Tensor[n_,__] -> Tensor[n,__] , + dot[Tensor[n_,__]] -> dot[Tensor[n,__]]}, 1]; + eqs = FilterRules[lookup[calc, Equations], patterns]; + calc3 = mapReplace[calc2, Equations, eqs]; + calc3 +]; + + + +(******************************************************************************) +(* Convert from ADMBase *) +(******************************************************************************) + +convertFromADMBaseCalc = +{ + Name -> thorn <> "_convertFromADMBase", + Schedule -> {"AT initial AFTER ADMBase_PostInitial"}, + ConditionalOnKeyword -> {"my_initial_data", "ADMBase"}, + Shorthands -> {g[la,lb], detg, gu[ua,ub], em4phi}, + Equations -> + { + g[la,lb] -> admg[la,lb], + detg -> detgExpr, + gu[ua,ub] -> 1/detg detgExpr MatrixInverse [g[ua,ub]], + + phi -> IfThen[conformalMethod==CMW, detg^(-1/6), Log[detg]/12], + em4phi -> IfThen[conformalMethod==CMW, 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], + + IfCCZ4[Theta -> 0] + } +}; + +convertFromADMBaseGammaCalc = +{ + Name -> thorn <> "_convertFromADMBaseGamma", + Schedule -> {"AT initial AFTER " <> thorn <> "_convertFromADMBase"}, + ConditionalOnKeyword -> {"my_initial_data", "ADMBase"}, + (* + Where -> InteriorNoSync, + *) + (* 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.) *) + Shorthands -> {dir[ua], + detgt, gtu[ua,ub], Gt[ua,lb,lc], theta}, + 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], + +(* + A -> - admdtalpha / (harmonicF alpha^harmonicN) (LapseAdvectionCoeff - 1), +*) + (* 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 Upwind[beta[ua], 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 Upwind[beta[ub], beta[ua], lb]), + 0] + } +}; + +(* Initialise the Gamma variables to 0. This is necessary with + multipatch because convertFromADMBaseGamma does not perform the + conversion in the boundary points, and the order in which symmetry + (interpatch) and outer boundary conditions is applied means that + points which are both interpatch and symmetry points are never + initialised. *) +initGammaCalc = +{ + Name -> thorn <> "_InitGamma", + Schedule -> {"AT initial BEFORE " <> thorn <> "_convertFromADMBaseGamma"}, + ConditionalOnKeyword -> {"my_initial_data", "ADMBase"}, + Where -> Everywhere, + Equations -> + { + Xt[ua] -> 0, + A -> 0, + B[ua] -> 0 + } +}; + + + +(******************************************************************************) +(* Convert to ADMBase *) +(******************************************************************************) + +convertToADMBaseCalc = +{ + Name -> thorn <> "_convertToADMBase", + Schedule -> {"IN " <> thorn <> "_convertToADMBaseGroup"}, + Where -> Everywhere, + Shorthands -> {e4phi}, + Equations -> + { + e4phi -> IfThen[conformalMethod==CMW, 1/phi^2, Exp[4 phi]], + admg[la,lb] -> e4phi gt[la,lb], + admK[la,lb] -> e4phi At[la,lb] + (1/3) admg[la,lb] trK, + admalpha -> alpha, + admbeta[ua] -> beta[ua] + } +}; + +convertToADMBaseDtLapseShiftCalc = +{ + Name -> thorn <> "_convertToADMBaseDtLapseShift", + Schedule -> {"IN " <> thorn <> "_convertToADMBaseGroup"}, + ConditionalOnKeyword -> {"dt_lapse_shift_method", "correct"}, + Where -> Interior, + Shorthands -> {dir[ua], detgt, gtu[ua,ub], eta, theta}, + Equations -> + { + dir[ua] -> Sign[beta[ua]], + + detgt -> 1 (* detgtExpr *), + (* This leads to simpler code... *) + gtu[ua,ub] -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]], + + eta -> etaExpr, + theta -> thetaExpr, + + (* see RHS *) +(* + admdtalpha -> - harmonicF alpha^harmonicN + ((1 - LapseAdvectionCoeff) A + LapseAdvectionCoeff trK) + + LapseAdvectionCoeff beta[ua] PDu[alpha,la], +*) + admdtalpha -> - harmonicF alpha^harmonicN + (+ LapseACoeff A + + ((1 - LapseACoeff) + (trK - IfCCZ4[2 Theta, 0]))) + + LapseAdvectionCoeff Upwind[beta[ua], alpha, la], + admdtbeta[ua] -> IfThen[harmonicShift, + - 1/2 gtu[ua,uj] phi alpha + (- 2 alpha PD[phi,lj] + + 2 phi PD[alpha,lj] + + gtu[uk,ul] phi alpha + (PD[gt[lk,ll],lj] - 2 PD[gt[lj,lk],ll])), + (* else *) + + theta ShiftGammaCoeff + (+ ShiftBCoeff B[ua] + + (1 - ShiftBCoeff) + (Xt[ua] - eta BetaDriver beta[ua]))] + + ShiftAdvectionCoeff Upwind[beta[ub], beta[ua], lb] + } +}; + +convertToADMBaseDtLapseShiftBoundaryCalc = +{ + Name -> thorn <> "_convertToADMBaseDtLapseShiftBoundary", + Schedule -> {"IN " <> thorn <> "_convertToADMBaseGroup"}, + ConditionalOnKeyword -> {"dt_lapse_shift_method", "correct"}, + Where -> BoundaryWithGhosts, + Shorthands -> {detgt, gtu[ua,ub], eta, theta}, + Equations -> + { + detgt -> 1 (* detgtExpr *), + (* This leads to simpler code... *) + gtu[ua,ub] -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]], + + eta -> etaExpr, + theta -> thetaExpr, + + (* see RHS, but omit derivatives near the boundary *) +(* + admdtalpha -> - harmonicF alpha^harmonicN + ((1 - LapseAdvectionCoeff) A + LapseAdvectionCoeff trK), +*) + admdtalpha -> - harmonicF alpha^harmonicN + (+ LapseACoeff A + + ((1 - LapseACoeff) + (trK - IfCCZ4[2 Theta, 0]))), + admdtbeta[ua] -> IfThen[harmonicShift, + 0, + (* else *) + + theta ShiftGammaCoeff + (+ ShiftBCoeff B[ua] + + (1 - ShiftBCoeff) + (Xt[ua] - eta BetaDriver beta[ua]))] + } +}; + +convertToADMBaseFakeDtLapseShiftCalc = +{ + Name -> thorn <> "_convertToADMBaseFakeDtLapseShift", + Schedule -> {"IN " <> thorn <> "_convertToADMBaseGroup"}, + ConditionalOnKeyword -> {"dt_lapse_shift_method", "noLapseShiftAdvection"}, + Where -> Everywhere, + Shorthands -> {detgt, gtu[ua,ub], eta, theta}, + Equations -> + { + detgt -> 1 (* detgtExpr *), + (* This leads to simpler code... *) + gtu[ua,ub] -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]], + + eta -> etaExpr, + theta -> thetaExpr, + + (* see RHS, but omit derivatives everywhere (which is wrong, but + faster, since it does not require synchronisation or boundary + conditions) *) +(* + admdtalpha -> - harmonicF alpha^harmonicN + ((1 - LapseAdvectionCoeff) A + LapseAdvectionCoeff trK), +*) + admdtalpha -> - harmonicF alpha^harmonicN + (+ LapseACoeff A + + ((1 - LapseACoeff) + (trK - IfCCZ4[2 Theta, 0]))), + admdtbeta[ua] -> IfThen[harmonicShift, + 0, + (* else *) + + theta ShiftGammaCoeff + (+ ShiftBCoeff B[ua] + + (1 - ShiftBCoeff) + (Xt[ua] - eta BetaDriver beta[ua]))] + } +}; + +(******************************************************************************) +(* Evolution equations *) +(******************************************************************************) + +evolCalc = +{ + Name -> thorn <> "_RHS", + Schedule -> {"IN " <> thorn <> "_evolCalcGroup"}, + (* + Where -> Interior, + *) + (* Synchronise the RHS grid functions after this routine, so that + the refinement boundaries are set correctly before applying the + radiative boundary conditions. *) + Where -> InteriorNoSync, + ConditionalOnKeyword -> {"RHS_split", "combined"}, + Shorthands -> {dir[ua], + detgt, gtu[ua,ub], + Gt[ua,lb,lc], Gtl[la,lb,lc], Gtlu[la,lb,uc], Xtn[ua], + Rt[la,lb], Rphi[la,lb], R[la,lb], + Atm[ua,lb], Atu[ua,ub], + e4phi, em4phi, cdphi[la], cdphi2[la,lb], g[la,lb], detg, + gu[ua,ub], Ats[la,lb], trAts, eta, theta, + rho, S[la], trS, fac1, fac2, dottrK, dotXt[ua], + epsdiss[ua], IfCCZ4[Z[ua]], IfCCZ4[dotTheta]}, + Equations -> + { + dir[ua] -> Sign[beta[ua]], + + detgt -> 1 (* detgtExpr *), + + (* This leads to simpler code... *) + gtu[ua,ub] -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]], + Gtl[la,lb,lc] -> 1/2 + (PD[gt[lb,la],lc] + PD[gt[lc,la],lb] - PD[gt[lb,lc],la]), + Gtlu[la,lb,uc] -> gtu[uc,ud] Gtl[la,lb,ld], + Gt[ua,lb,lc] -> gtu[ua,ud] Gtl[ld,lb,lc], + + (* The conformal connection functions calculated from the conformal metric, + used instead of Xt where no derivatives of Xt are taken *) + Xtn[ui] -> gtu[uj,uk] Gt[ui,lj,lk], + + e4phi -> IfThen[conformalMethod==CMW, 1/phi^2, Exp[4 phi]], + em4phi -> 1 / e4phi, + g[la,lb] -> e4phi gt[la,lb], + detg -> detgExpr, + gu[ua,ub] -> em4phi gtu[ua,ub], + + (* The Z quantities *) + (* gr-qc:1106.2254 (2011), eqn. (23) *) + IfCCZ4[ + Z[ud] -> (1/2) gu[ua,ud] (- PD[gt[la,lb],lc] gtu[ub,uc] + gt[la,lc] Xt[uc]) + ], + + (* PRD 62, 044034 (2000), eqn. (18) *) + (* Adding Z term by changing Xtn to Xt *) + Rt[li,lj] -> - (1/2) gtu[ul,um] PD[gt[li,lj],ll,lm] + + (1/2) gt[lk,li] PD[Xt[uk],lj] + + (1/2) gt[lk,lj] PD[Xt[uk],li] + + (1/2) Xtn[uk] Gtl[li,lj,lk] + + (1/2) Xtn[uk] Gtl[lj,li,lk] + + (+ Gt[uk,li,ll] Gtlu[lj,lk,ul] + + Gt[uk,lj,ll] Gtlu[li,lk,ul] + + Gt[uk,li,ll] Gtlu[lk,lj,ul]), + + fac1 -> IfThen[conformalMethod==CMW, -1/(2 phi), 1], + cdphi[la] -> fac1 CDt[phi,la], + fac2 -> IfThen[conformalMethod==CMW, 1/(2 phi^2), 0], + cdphi2[la,lb] -> fac1 CDt[phi,la,lb] + fac2 CDt[phi,la] CDt[phi,lb], + + (* PRD 62, 044034 (2000), eqn. (15) *) + Rphi[li,lj] -> - 2 cdphi2[lj,li] + - 2 gt[li,lj] gtu[ul,un] cdphi2[ll,ln] + + 4 cdphi[li] cdphi[lj] + - 4 gt[li,lj] gtu[ul,un] cdphi[ln] cdphi[ll], + + Atm[ua,lb] -> gtu[ua,uc] At[lc,lb], + Atu[ua,ub] -> Atm[ua,lc] gtu[ub,uc], + + R[la,lb] -> Rt[la,lb] + Rphi[la,lb], + IfCCZ4[ + R[la,lb] -> R[la,lb] + (2/phi) (+ g[la,lc] Z[uc] PD[phi,lb] + + g[lb,lc] Z[uc] PD[phi,la] - g[la,lb] Z[uc] PD[phi,lc]) + + e4phi Z[uc] PD[gt[la,lb],lc] + ], + + (* Matter terms *) + + (* rho = n^a n^b T_ab *) + rho -> addMatter + (1/alpha^2 (T00 - 2 beta[ui] T0[li] + beta[ui] beta[uj] T[li,lj])), + + (* S_i = -p^a_i n^b T_ab, where p^a_i = delta^a_i + n^a n_i *) + S[li] -> addMatter (-1/alpha (T0[li] - beta[uj] T[li,lj])), + + (* trS = gamma^ij T_ij *) + trS -> addMatter (em4phi gtu[ui,uj] T[li,lj]), + + (* RHS terms *) + + (* PRD 62, 044034 (2000), eqn. (10) *) + (* PRD 67 084023 (2003), eqn. (16) and (23) *) + dot[phi] -> IfThen[conformalMethod==CMW, 1/3 phi, -1/6] + (alpha trK - PD[beta[ua],la]), + + (* PRD 62, 044034 (2000), eqn. (9) *) + (* gr-qc:1106.2254 (2011), eqn. (14) *) + (* removing trA from Aij ensures that detg = 1 *) + dot[gt[la,lb]] -> - 2 alpha (At[la,lb] - IfCCZ4[(1/3) At[lc,ld] gtu[uc,ud] gt[la,lb], 0]) + + gt[la,lc] PD[beta[uc],lb] + gt[lb,lc] PD[beta[uc],la] + - (2/3) gt[la,lb] PD[beta[uc],lc], + (* PRD 62, 044034 (2000), eqn. (20) *) + (* PRD 67 084023 (2003), eqn (26) *) + (* gr-qc:1106.2254 (2011), eqn. (19) *) + (* Adding Z terms by changing Xtn to Xt, + also adding extra Z and Theta terms *) + dotXt[ui] -> - 2 Atu[ui,uj] PD[alpha,lj] + + 2 alpha (+ Gt[ui,lj,lk] Atu[uk,uj] + - (2/3) gtu[ui,uj] PD[trK,lj] + + 6 Atu[ui,uj] cdphi[lj]) + + gtu[uj,ul] PD[beta[ui],lj,ll] + + (1/3) gtu[ui,uj] PD[beta[ul],lj,ll] + - Xtn[uj] PD[beta[ui],lj] + + (2/3) Xtn[ui] PD[beta[uj],lj] + + IfCCZ4[ + + GammaShift 2 e4phi (- Z[uj] PD[beta[ui],lj] + + (2/3) Z[ui] PD[beta[uj],lj]) + - (4/3) alpha e4phi Z[ui] trK + + 2 gtu[ui,uj] (+ alpha PD[Theta,lj] + - Theta PD[alpha,lj]) + - 2 alpha e4phi dampk1 Z[ui], + 0] + (* Equation (4.28) in Baumgarte & Shapiro (Phys. Rept. 376 (2003) 41-131) *) + + addMatter (- 16 Pi alpha gtu[ui,uj] S[lj]), + dot[Xt[ui]] -> dotXt[ui], + + (* gr-qc:1106.2254 (2011), eqn. (18) *) + IfCCZ4[ + dotTheta -> + - PD[alpha,la] Z[ua] - dampk1 (2 + dampk2) alpha Theta + + (1/2) alpha (gu[ua,ub] R[la,lb] - Atm[ua,lb] Atm[ub,la] + (2/3) trK^2 - 2 trK Theta) + + addMatter (- 8 Pi alpha rho) + ], + + IfCCZ4[ + dot[Theta] -> dotTheta + ], + + (* PRD 62, 044034 (2000), eqn. (11) *) + (* gr-qc:1106.2254 (2011), eqn. (17) *) + (* Adding the RHS of Theta to K, because K_Z4 = K_BSSN + 2 Theta *) + (* Also adding the Z term, as it has to cancel with the one in Theta *) + dottrK -> - em4phi ( gtu[ua,ub] ( PD[alpha,la,lb] + + 2 cdphi[la] PD[alpha,lb] ) + - Xtn[ua] PD[alpha,la] ) + + alpha (Atm[ua,lb] Atm[ub,la] + (1/3) trK^2) + + IfCCZ4[ + + 2 dotTheta + 2 PD[alpha,la] Z[ua] + + dampk1 (1 - dampk2) alpha Theta, + 0] + (* Equation (4.21) in Baumgarte & Shapiro (Phys. Rept. 376 (2003) 41-131) *) + + addMatter (4 Pi alpha (rho + trS)), + dot[trK] -> dottrK, + + (* PRD 62, 044034 (2000), eqn. (12) *) + (* TODO: Should we use the Hamiltonian constraint to make Rij tracefree? *) + (* gr-qc:1106.2254 (2011), eqn. (15) *) + (* Adding Z terms in the Ricci and Theta terms *) + Ats[la,lb] -> - CDt[alpha,la,lb] + + + 2 (PD[alpha,la] cdphi[lb] + PD[alpha,lb] cdphi[la] ) + + alpha R[la,lb], + trAts -> gu[ua,ub] Ats[la,lb], + dot[At[la,lb]] -> + em4phi (+ Ats[la,lb] - (1/3) g[la,lb] trAts ) + + alpha (+ ((trK - IfCCZ4[2 Theta, 0]) + At[la,lb]) + - 2 At[la,lc] Atm[uc,lb]) + + At[la,lc] PD[beta[uc],lb] + At[lb,lc] PD[beta[uc],la] + - (2/3) At[la,lb] PD[beta[uc],lc] + (* Equation (4.23) in Baumgarte & Shapiro (Phys. Rept. 376 (2003) 41-131) *) + + addMatter (- em4phi alpha 8 Pi + (T[la,lb] - (1/3) g[la,lb] trS)), + + (* dot[alpha] -> - harmonicF alpha^harmonicN trK, *) + (* dot[alpha] -> - harmonicF alpha^harmonicN A + Lie[alpha, beta], *) +(* + dot[alpha] -> - harmonicF alpha^harmonicN ( + (1 - LapseAdvectionCoeff) A + LapseAdvectionCoeff trK) + + LapseAdvectionCoeff beta[ua] PDu[alpha,la], + dot[A] -> (1 - LapseAdvectionCoeff) (dottrK - AlphaDriver A), +*) + dot[alpha] -> - harmonicF alpha^harmonicN + (+ LapseACoeff A + + ((1 - LapseACoeff) + (+ trK - IfCCZ4[2 Theta, 0] + + AlphaDriver (alpha - 1)))), + + dot[A] -> + (LapseACoeff + (+ dottrK - IfCCZ4[2 dotTheta, 0] + - AlphaDriver A)), + + eta -> etaExpr, + theta -> thetaExpr, + + (* dot[beta[ua]] -> eta Xt[ua], *) + (* dot[beta[ua]] -> ShiftGammaCoeff alpha^ShiftAlphaPower B[ua], *) + dot[beta[ua]] -> IfThen[harmonicShift, + - 1/2 gtu[ua,uj] phi alpha + (- 2 alpha PD[phi,lj] + + 2 phi PD[alpha,lj] + + gtu[uk,ul] phi alpha + (PD[gt[lk,ll],lj] - 2 PD[gt[lj,lk],ll])), + (* else *) + + theta ShiftGammaCoeff + (+ ShiftBCoeff B[ua] + + (1 - ShiftBCoeff) + (Xt[ua] - eta BetaDriver beta[ua]))], + + dot[B[ua]] -> + ShiftBCoeff (dotXt[ua] - eta BetaDriver B[ua]) + (* Note that this dotXt[ua] is not yet \partial_t \tilde \Gamma^i, because the + advection term has not yet been added. It is actually + \partial_t \tilde \Gamma^i - \beta^j \partial_j \tilde \Gamma^i *) + } +}; + +advectCalc = +{ + Name -> thorn <> "_Advect", + Schedule -> {"IN " <> thorn <> "_evolCalcGroup " <> + "AFTER (" <> thorn <> "_RHS " <> thorn <> "_RHS1 " <> thorn <> "_RHS2)"}, + (* + Where -> Interior, + *) + (* Synchronise the RHS grid functions after this routine, so that + the refinement boundaries are set correctly before applying the + radiative boundary conditions. *) + Where -> InteriorNoSync, + ConditionalOnKeyword -> {"advection_split", "combined"}, + Shorthands -> {dir[ua]}, + Equations -> + { + dir[ua] -> Sign[beta[ua]], + + dot[phi] -> dot[phi] + Upwind[beta[ua], phi, la], + + dot[gt[la,lb]] -> dot[gt[la,lb]] + Upwind[beta[uc], gt[la,lb], lc], + + dot[Xt[ui]] -> dot[Xt[ui]] + Upwind[beta[uj], Xt[ui], lj], + + IfCCZ4[ + dot[Theta] -> dot[Theta] + Upwind[beta[ua], Theta, la] + ], + + dot[trK] -> dot[trK] + Upwind[beta[ua], trK, la], + + dot[At[la,lb]] -> dot[At[la,lb]] + Upwind[beta[uc], At[la,lb], lc], + + dot[alpha] -> dot[alpha] + + LapseAdvectionCoeff Upwind[beta[ua], alpha, la], + + dot[A] -> dot[A] + + LapseACoeff ( + + LapseAdvectionCoeff Upwind[beta[ua], A, la] + + (1 - LapseAdvectionCoeff) Upwind[beta[ua], trK, la]), + + dot[beta[ua]] -> dot[beta[ua]] + + ShiftAdvectionCoeff Upwind[beta[ub], beta[ua], lb], + + dot[B[ua]] -> dot[B[ua]] + + ShiftBCoeff ( + + ShiftAdvectionCoeff Upwind[beta[ub], B[ua], lb] + + ((1 - ShiftAdvectionCoeff) + Upwind[beta[ub], Xt[ua], lb])) + (* Note that the advection term \beta^j \partial_j \tilde \Gamma^i is not + subtracted here when ShiftAdvectionCoefficient == 1 because it was + implicitly subtracted before (see comment in previous calculation of + dot[B[ua]]. *) + } +}; + +varsNames = { + {"phi", dot[phi]}, + {"gt", dot[gt[la,lb]]}, + {"Xt", dot[Xt[ui]]}, + {"trK", dot[trK]}, + {"At", dot[At[la,lb]]}, + {"alpha", dot[alpha]}, + {"A", dot[A]}, + {"beta", dot[beta[ua]]}, + {"B", dot[B[ua]]}, + IfCCZ4[{"Theta", dot[Theta]}] + }; + +advectCalcs = Map[ + PartialCalculation[advectCalc, "_"<>ToString[First[#]], + {ConditionalOnKeyword -> {"advection_split", + "per variable"}}, + {Last[#]}]&, + varsNames]; + +evolCalc1 = PartialCalculation[evolCalc, "1", + { + ConditionalOnKeyword -> {"RHS_split", "split At"} + }, + { + dot[phi], + dot[gt[la,lb]], + dot[Xt[ui]], + dot[trK], + dot[alpha], + dot[A], + dot[beta[ua]], + dot[B[ua]], + IfCCZ4[dot[Theta]] + }]; + +evolCalc2 = PartialCalculation[evolCalc, "2", + { + ConditionalOnKeyword -> {"RHS_split", "split At"} + }, + { + dot[At[la,lb]] + }]; + +dissCalc = +{ + Name -> thorn <> "_Dissipation", + Schedule -> {"IN " <> thorn <> "_evolCalcGroup " <> + "AFTER (" <> thorn <> "_RHS1 " <> thorn <> "_RHS2)"}, + ConditionalOnKeyword -> {"apply_dissipation", "always"}, + Where -> InteriorNoSync, + Shorthands -> {epsdiss[ua]}, + Equations -> + { + epsdiss[ua] -> EpsDiss, + Sequence@@Table[ + dot[var] -> dot[var] + epsdiss[ux] PDdiss[var,lx], + {var, {phi, gt[la,lb], Xt[ui], IfCCZ4[Theta], trK, At[la,lb], + alpha, A, beta[ua], B[ua]}}] + } +}; + +dissCalcs = +Table[ +{ + Name -> thorn <> "_Dissipation_" <> ToString[var /. {Tensor[n_,__] -> n}], + Schedule -> {"IN " <> thorn <> "_evolCalcGroup " <> + "AFTER (" <> thorn <> "_RHS1 " <> thorn <> "_RHS2)"}, + ConditionalOnKeyword -> {"apply_dissipation", "always"}, + Where -> InteriorNoSync, + Shorthands -> {epsdiss[ua]}, + Equations -> + { + epsdiss[ua] -> EpsDiss, + dot[var] -> dot[var] + epsdiss[ux] PDdiss[var,lx] + } +}, + {var, {phi, gt[la,lb], Xt[ui], IfCCZ4[Theta], trK, At[la,lb], + alpha, A, beta[ua], B[ua]}} +]; + +RHSStaticBoundaryCalc = +{ + Name -> thorn <> "_RHSStaticBoundary", + Schedule -> {"IN MoL_CalcRHS"}, + ConditionalOnKeyword -> {"my_rhs_boundary_condition", "static"}, + Where -> Boundary, + Equations -> + { + dot[phi] -> 0, + dot[gt[la,lb]] -> 0, + dot[trK] -> 0, + dot[At[la,lb]] -> 0, + dot[Xt[ua]] -> 0, + dot[alpha] -> 0, + dot[A] -> 0, + dot[beta[ua]] -> 0, + dot[B[ua]] -> 0, + IfCCZ4[dot[Theta] -> 0] + } +}; + +(* Initialise the RHS variables in analysis in case they are going to + be output - the noninterior points cannot be filled, so we define + them to be zero *) +initRHSCalc = +{ + Name -> thorn <> "_InitRHS", + Schedule -> {"AT analysis BEFORE " <> thorn <> "_evolCalcGroup"}, + Where -> Everywhere, + Equations -> + { + dot[phi] -> 0, + dot[gt[la,lb]] -> 0, + dot[trK] -> 0, + dot[At[la,lb]] -> 0, + dot[Xt[ua]] -> 0, + dot[alpha] -> 0, + dot[A] -> 0, + dot[beta[ua]] -> 0, + dot[B[ua]] -> 0, + IfCCZ4[dot[Theta] -> 0] + } +}; + +RHSRadiativeBoundaryCalc = +{ + Name -> thorn <> "_RHSRadiativeBoundary", + Schedule -> {"IN MoL_CalcRHS"}, + ConditionalOnKeyword -> {"my_rhs_boundary_condition", "radiative"}, + Where -> Boundary, + Shorthands -> {dir[ua], + detgt, gtu[ua,ub], em4phi, gu[ua,ub], + nn[la], nu[ua], nlen, nlen2, su[ua], + vg}, + Equations -> + { + dir[ua] -> Sign[normal[ua]], + + detgt -> 1 (* detgtExpr *), + gtu[ua,ub] -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]], + em4phi -> IfThen[conformalMethod==CMW, phi^2, Exp[-4 phi]], + gu[ua,ub] -> em4phi gtu[ua,ub], + + nn[la] -> Euc[la,lb] normal[ub], + nu[ua] -> gu[ua,ub] nn[lb], + nlen2 -> nu[ua] nn[la], + nlen -> Sqrt[nlen2], + su[ua] -> nu[ua] / nlen, + + vg -> Sqrt[harmonicF], + + dot[phi] -> - vg su[uc] PDo[phi ,lc], + dot[gt[la,lb]] -> - su[uc] PDo[gt[la,lb],lc], + dot[trK] -> - vg su[uc] PDo[trK ,lc], + dot[At[la,lb]] -> - su[uc] PDo[At[la,lb],lc], + dot[Xt[ua]] -> - su[uc] PDo[Xt[ua] ,lc], + dot[alpha] -> - vg su[uc] PDo[alpha ,lc], + dot[A] -> - vg su[uc] PDo[A ,lc], + dot[beta[ua]] -> - su[uc] PDo[beta[ua] ,lc], + dot[B[ua]] -> - su[uc] PDo[B[ua] ,lc], + IfCCZ4[ + dot[Theta] -> - vg su[uc] PDo[Theta ,lc] + ] + } +}; + +enforceCalc = +{ + Name -> thorn <> "_enforce", + Schedule -> {"IN MoL_PostStepModify"}, + Shorthands -> {detgt, gtu[ua,ub], trAt}, + Equations -> + { + (* The following comment is still interesting, but is not correct + any more since it is now scheduled in MoL_PostStepModify instead: + + Enforcing the constraints needs to be a projection, because it + is applied in MoL_PostStep and may thus be applied multiple + times, not only during time evolution. Therefore detgt has to + be calculated correctly, without assuming that det gt_ij = 1, + which is not always the case (since we don't enforce it). On + the other hand, this may not be so important... *) + detgt -> 1 (* detgtExpr *), + gtu[ua,ub] -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]], + + trAt -> gtu[ua,ub] At[la,lb], + + At[la,lb] -> At[la,lb] - (1/3) gt[la,lb] trAt, + + alpha -> Max[alpha, MinimumLapse] + } +}; + +(******************************************************************************) +(* Boundary conditions *) +(******************************************************************************) + +boundaryCalc = +{ + Name -> thorn <> "_boundary", + Schedule -> {"IN MoL_PostStep"}, + ConditionalOnKeyword -> {"my_boundary_condition", "Minkowski"}, + Where -> BoundaryWithGhosts, + Equations -> + { + phi -> IfThen[conformalMethod==CMW, 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, + IfCCZ4[Theta -> 0] + } +}; + +(******************************************************************************) +(* Constraint equations *) +(******************************************************************************) + +constraintsCalc = +{ + Name -> thorn <> "_constraints", + Schedule -> Automatic, + After -> "MoL_PostStep", + Where -> Interior, + Shorthands -> {detgt, ddetgt[la], gtu[ua,ub], Z[ua], + Gt[ua,lb,lc], Gtl[la,lb,lc], Gtlu[la,lb,uc], Xtn[ua], + e4phi, em4phi, + g[la,lb], detg, gu[ua,ub], ddetg[la], G[ua,lb,lc], + Rt[la,lb], Rphi[la,lb], R[la,lb], trR, Atm[ua,lb], + gK[la,lb,lc], cdphi[la], cdphi2[la,lb], + rho, S[la], fac1, fac2}, + Equations -> + { + detgt -> 1 (* detgtExpr *), + ddetgt[la] -> 0 (* ddetgtExpr[la] *), + + (* This leads to simpler code... *) + gtu[ua,ub] -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]], + Gtl[la,lb,lc] -> 1/2 + (PD[gt[lb,la],lc] + PD[gt[lc,la],lb] - PD[gt[lb,lc],la]), + Gtlu[la,lb,uc] -> gtu[uc,ud] Gtl[la,lb,ld], + Gt[ua,lb,lc] -> gtu[ua,ud] Gtl[ld,lb,lc], + + (* The conformal connection functions calculated from the conformal metric, + used instead of Xt where no derivatives of Xt are taken *) + Xtn[ui] -> gtu[uj,uk] Gt[ui,lj,lk], + + e4phi -> IfThen[conformalMethod==CMW, 1/phi^2, Exp[4 phi]], + em4phi -> 1 / e4phi, + g[la,lb] -> e4phi gt[la,lb], + detg -> e4phi^3, + gu[ua,ub] -> em4phi gtu[ua,ub], + + (* The Z quantities *) + IfCCZ4[ + Z[ud] -> (1/2) gu[ua,ud] (- PD[gt[la,lb],lc] gtu[ub,uc] + gt[la,lc] Xt[uc]) + ], + + (* PRD 62, 044034 (2000), eqn. (18) *) + Rt[li,lj] -> - (1/2) gtu[ul,um] PD[gt[li,lj],ll,lm] + + (1/2) gt[lk,li] PD[Xt[uk],lj] + + (1/2) gt[lk,lj] PD[Xt[uk],li] + + (1/2) Xtn[uk] Gtl[li,lj,lk] + + (1/2) Xtn[uk] Gtl[lj,li,lk] + + (+ Gt[uk,li,ll] Gtlu[lj,lk,ul] + + Gt[uk,lj,ll] Gtlu[li,lk,ul] + + Gt[uk,li,ll] Gtlu[lk,lj,ul]), + + (* From the long turducken paper. + This expression seems to give the same result as the one from 044034. *) + (* TODO: symmetrise correctly: (ij) = (1/2) [i+j] *) +(* + Rt[li,lj] -> - (1/2) gtu[uk,ul] PD[gt[li,lj],lk,ll] + + gt[lk,li] PD[Xt[uk],lj] + gt[lk,lj] PD[Xt[uk],li] + + gt[li,ln] Gt[un,lj,lk] gtu[um,ua] gtu[uk,ub] PD[gt[la,lb],lm] + + gt[lj,ln] Gt[un,li,lk] gtu[um,ua] gtu[uk,ub] PD[gt[la,lb],lm] + + gtu[ul,us] (+ 2 Gt[uk,ll,li] gt[lj,ln] Gt[un,lk,ls] + + 2 Gt[uk,ll,lj] gt[li,ln] Gt[un,lk,ls] + + Gt[uk,li,ls] gt[lk,ln] Gt[un,ll,lj]), +*) + + (* Below would be a straightforward calculation, + without taking any Gamma^i into account. + This expression gives a different answer! *) +(* + Rt[la,lb] -> + Gt[u1,l2,la] Gt[l1,lb,u2] - Gt[u1,la,lb] Gt[l1,l2,u2] + + 1/2 gtu[u1,u2] (- PD[gt[l1,l2],la,lb] + PD[gt[l1,la],l2,lb] + - PD[gt[la,lb],l1,l2] + PD[gt[l2,lb],l1,la]), +*) + + fac1 -> IfThen[conformalMethod==CMW, -1/(2 phi), 1], + cdphi[la] -> fac1 CDt[phi,la], + fac2 -> IfThen[conformalMethod==CMW, 1/(2 phi^2), 0], + cdphi2[la,lb] -> fac1 CDt[phi,la,lb] + fac2 CDt[phi,la] CDt[phi,lb], + + (* PRD 62, 044034 (2000), eqn. (15) *) + Rphi[li,lj] -> - 2 cdphi2[lj,li] + - 2 gt[li,lj] gtu[ul,un] cdphi2[ll,ln] + + 4 cdphi[li] cdphi[lj] + - 4 gt[li,lj] gtu[ul,un] cdphi[ln] cdphi[ll], + + (* ddetg[la] -> PD[e4phi detg,la], *) + ddetg[la] -> e4phi ddetgt[la] + 4 detgt e4phi PD[phi,la], + (* TODO: check this equation, maybe simplify it by omitting ddetg *) + G[ua,lb,lc] -> Gt[ua,lb,lc] + + 1/(2 detg) (+ KD[ua,lb] ddetg[lc] + KD[ua,lc] ddetg[lb] + - (1/3) g[lb,lc] gu[ua,ud] ddetg[ld]), + + R[la,lb] -> + Rt[la,lb] + Rphi[la,lb], + + IfCCZ4[ + R[la,lb] -> R[la, lb] + (2/phi) (+ g[la,lc] Z[uc] PD[phi,lb] + + g[lb,lc] Z[uc] PD[phi,la] - g[la,lb] Z[uc] PD[phi,lc]) + + e4phi Z[uc] PD[gt[la,lb],lc] + ], + + trR -> gu[ua,ub] R[la,lb], + + (* K[la,lb] -> e4phi At[la,lb] + (1/3) g[la,lb] trK, *) + (* Km[ua,lb] -> gu[ua,uc] K[lc,lb], *) + Atm[ua,lb] -> gtu[ua,uc] At[lc,lb], + + (* Matter terms *) + + (* rho = n^a n^b T_ab *) + rho -> 1/alpha^2 (T00 - 2 beta[ui] T0[li] + beta[ui] beta[uj] T[li,lj]), + + (* S_i = -p^a_i n^b T_ab, where p^a_i = delta^a_i + n^a n_i *) + S[li] -> -1/alpha (T0[li] - beta[uj] T[li,lj]), + + (* Constraints *) + + (* H -> trR - Km[ua,lb] Km[ub,la] + trK^2, *) + (* PRD 67, 084023 (2003), eqn. (19) *) + H -> trR - Atm[ua,lb] Atm[ub,la] + (2/3) trK^2 - addMatter 16 Pi rho, + + (* gK[la,lb,lc] -> CD[K[la,lb],lc], *) +(* gK[la,lb,lc] -> + 4 e4phi PD[phi,lc] At[la,lb] + e4phi CD[At[la,lb],lc] + + (1/3) g[la,lb] PD[trK,lc], + M[la] -> gu[ub,uc] (gK[lc,la,lb] - gK[lc,lb,la]), *) + + M[li] -> + gtu[uj,uk] (CDt[At[li,lj],lk] + 6 At[li,lj] cdphi[lk]) + - (2/3) PD[trK,li] + - addMatter 8 Pi S[li], + (* TODO: use PRD 67, 084023 (2003), eqn. (20) *) + + (* det gamma-tilde *) + cS -> Log[detgt], + + (* Gamma constraint *) + cXt[ua] -> gtu[ub,uc] Gt[ua,lb,lc] - Xt[ua], + + (* trace A-tilde *) + cA -> gtu[ua,ub] At[la,lb] + } +}; + +constraintsCalc1 = PartialCalculation[constraintsCalc, "1", + {}, + { + H + }]; + +constraintsCalc2 = PartialCalculation[constraintsCalc, "2", + {}, + { + M[li], + cS, + cXt[ua], + cA + }]; + +(******************************************************************************) +(* Implementations *) +(******************************************************************************) + +inheritedImplementations = + Join[{"ADMBase"}, + If [addMatter!=0, {"TmunuBase"}, {}]]; + +(******************************************************************************) +(* Parameters *) +(******************************************************************************) + +inheritedKeywordParameters = {}; + +extendedKeywordParameters = +{ + { + Name -> "ADMBase::evolution_method", + AllowedValues -> {thorn} + }, + { + Name -> "ADMBase::lapse_evolution_method", + AllowedValues -> {thorn} + }, + { + Name -> "ADMBase::shift_evolution_method", + AllowedValues -> {thorn} + }, + { + Name -> "ADMBase::dtlapse_evolution_method", + AllowedValues -> {thorn} + }, + { + Name -> "ADMBase::dtshift_evolution_method", + AllowedValues -> {thorn} + } +}; + +keywordParameters = +{ + { + Name -> "my_initial_data", + (* Visibility -> "restricted", *) + (* Description -> "ddd", *) + AllowedValues -> {"ADMBase", "Minkowski"}, + Default -> "ADMBase" + }, + { + Name -> "my_initial_boundary_condition", + Visibility -> "restricted", + (* Description -> "ddd", *) + AllowedValues -> {"none"}, + Default -> "none" + }, + { + Name -> "my_rhs_boundary_condition", + Visibility -> "restricted", + (* Description -> "ddd", *) + AllowedValues -> {"none", "static", "radiative"}, + Default -> "none" + }, + { + Name -> "my_boundary_condition", + (* Visibility -> "restricted", *) + (* Description -> "ddd", *) + AllowedValues -> {"none", "Minkowski"}, + Default -> "none" + }, + { + Name -> "calculate_ADMBase_variables_at", + Visibility -> "restricted", + (* Description -> "ddd", *) + AllowedValues -> {"MoL_PostStep", "CCTK_EVOL", "CCTK_ANALYSIS"}, + Default -> "MoL_PostStep" + }, + { + Name -> "UseSpatialBetaDriver_UNUSED", + Visibility -> "restricted", + (* Description -> "ddd", *) + AllowedValues -> {"no", "yes"}, + Default -> "no" + }, + { + Name -> "dt_lapse_shift_method", + Description -> "Treatment of ADMBase dtlapse and dtshift", + AllowedValues -> {"correct", + "noLapseShiftAdvection" (* omit lapse and shift advection terms (faster) *) + }, + Default -> "correct" + }, + { + Name -> "apply_dissipation", + Description -> "Whether to apply dissipation to the RHSs", + AllowedValues -> {"always", + "never" (* yes and no keyword values confuse Cactus, and Kranc + doesn't support boolean parameters *) + }, + Default -> "never" + }, + { + Name -> "RHS_split", + Description -> "How to split RHS calculation", + AllowedValues -> {"combined", + "split At"}, + Default -> "split At" + }, + { + Name -> "advection_split", + Description -> "How to split advection calculation", + AllowedValues -> {"combined", + "per variable"}, + Default -> "combined" + } +}; + +intParameters = +{ + { + Name -> harmonicN, + Description -> "d/dt alpha = - f alpha^n K (harmonic=2, 1+log=1)", + Default -> 2 + }, + { + Name -> ShiftAlphaPower, + Default -> 0 + }, + { + Name -> conformalMethod, + Description -> "Treatment of conformal factor", + AllowedValues -> {{Value -> "0", Description -> "phi method"}, + {Value -> "1", Description -> "W method"}}, + Default -> 0 + }, + { + Name -> fdOrder, + Default -> derivOrder, + AllowedValues -> {2,4,6,8} + }, + { + Name -> harmonicShift, + Description -> "Whether to use the harmonic shift", + AllowedValues -> {{Value -> "0", Description -> "Gamma driver shift"}, + {Value -> "1", Description -> "Harmonic shift"}}, + Default -> 0 + } +}; + +realParameters = +{ + IfCCZ4[{ + Name -> GammaShift, + Description -> "Covariant shift term in Gamma", + Default -> 0.5 + }], + IfCCZ4[{ + Name -> dampk1, + Description -> "CCZ4 damping term 1 for Theta and Z", + Default -> 0 + }], + IfCCZ4[{ + Name -> dampk2, + Description -> "CCZ4 damping term 2 for Theta and Z", + Default -> 0 + }], + { + Name -> LapseACoeff, + Description -> "Whether to evolve A in time", + Default -> 0 + }, + { + Name -> harmonicF, + Description -> "d/dt alpha = - f alpha^n K (harmonic=1, 1+log=2)", + Default -> 1 + }, + { + Name -> AlphaDriver, + Default -> 0 + }, + { + Name -> ShiftBCoeff, + Description -> "Whether to evolve B^i in time", + Default -> 1 + }, + { + Name -> ShiftGammaCoeff, + Default -> 0 + }, + { + Name -> BetaDriver, + Default -> 0 + }, + { + Name -> LapseAdvectionCoeff, + Description -> "Factor in front of the lapse advection terms in 1+log", + Default -> 1 + }, + { + Name -> ShiftAdvectionCoeff, + Description -> "Factor in front of the shift advection terms in gamma driver", + Default -> 1 + }, + { + Name -> MinimumLapse, + Description -> "Minimum value of the lapse function", + Default -> -1 + }, + { + Name -> SpatialBetaDriverRadius, + Description -> "Radius at which the BetaDriver starts to be reduced", + AllowedValues -> {{Value -> "(0:*", Description -> "Positive"}}, + Default -> 10^12 + }, + { + Name -> SpatialShiftGammaCoeffRadius, + Description -> "Radius at which the ShiftGammaCoefficient starts to be reduced", + AllowedValues -> {{Value -> "(0:*", Description -> "Positive"}}, + Default -> 10^12 + }, + { + Name -> EpsDiss, + Description -> "Dissipation strength", + AllowedValues -> {{Value -> "(0:*", Description -> "Positive"}}, + Default -> 0 + } +}; + +(******************************************************************************) +(* Construct the thorns *) +(******************************************************************************) + +calculations = +Join[ +{ + initialCalc, + convertFromADMBaseCalc, + initGammaCalc, + convertFromADMBaseGammaCalc, + evolCalc, + evolCalc1, evolCalc2, + dissCalc, + advectCalc, + (*advectCalcs,*) + initRHSCalc, + RHSStaticBoundaryCalc, + (* RHSRadiativeBoundaryCalc, *) + enforceCalc, + boundaryCalc, + convertToADMBaseCalc, + convertToADMBaseDtLapseShiftCalc, + convertToADMBaseDtLapseShiftBoundaryCalc, + convertToADMBaseFakeDtLapseShiftCalc, + (* constraintsCalc, *) + constraintsCalc1, constraintsCalc2 +}, + advectCalcs + (*dissCalcs*) +]; + +CreateKrancThornTT [groups, ".", thorn, + Calculations -> calculations, + DeclaredGroups -> declaredGroupNames, + PartialDerivatives -> derivatives, + EvolutionTimelevels -> evolutionTimelevels, + DefaultEvolutionTimelevels -> 3, + UseJacobian -> True, + UseLoopControl -> True, + UseVectors -> useVectors, + UseOpenCL -> useOpenCL, + InheritedImplementations -> inheritedImplementations, + InheritedKeywordParameters -> inheritedKeywordParameters, + ExtendedKeywordParameters -> extendedKeywordParameters, + KeywordParameters -> keywordParameters, + IntParameters -> intParameters, + RealParameters -> realParameters +]; + +]; + + + +(******************************************************************************) +(* Options *) +(******************************************************************************) + +(* These are the arguments to createComment: + - derivative order: 2, 4, 6, 8, ... + - useJacobian: False or True + - split upwind derivatives: False or True + - use vectorisation: False or True + - use OpenCL: 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) + - thorn base name +*) + +createCode[4, False, True, True , False, 3, 1, "BSSN"]; +createCode[4, False, True, False, False, 3, 1, "BSSN"]; +createCode[4, False, True, True , True , 3, 1, "BSSN"]; + +createCode[4, False, True, True , False, 3, 1, "CCZ4"]; + From e71809861e3579a680c011cac3a5eb45006149e9 Mon Sep 17 00:00:00 2001 From: Alex Johnson Date: Wed, 28 Jun 2017 17:38:48 +0530 Subject: [PATCH 2/3] OTWO-3061 Octave comments nested within mathematica When octave comments were nested within mathematica comments, the parser incorrectly detected the file as Octave. To fix this we moved around an If/Else block which placed mathematica at the bottom of the conditional. We have also tried adding mathematica comments within Octave comments to see if the reverse was true. We found that Octave parsing worked correctly without any further change in code. --- src/detector.c | 8 +++----- test/expected_dir/mathematica1.m | 4 ++-- test/expected_dir/octave1.m | 6 +++--- test/src_dir/mathematica1.m | 4 ++-- test/src_dir/octave1.m | 6 +++--- 5 files changed, 13 insertions(+), 15 deletions(-) diff --git a/src/detector.c b/src/detector.c index 414fea4..fa09f49 100644 --- a/src/detector.c +++ b/src/detector.c @@ -788,12 +788,10 @@ const char *disambiguate_m(SourceFile *sourcefile) { else if (objective_c_score > matlab_score && objective_c_score > mathematica_score) return LANG_OBJECTIVE_C; - else if (octave_syntax_detected) - return LANG_OCTAVE; - else if (matlab_score > mathematica_score) - return LANG_MATLAB; - else + else if (mathematica_score > matlab_score) return LANG_MATHEMATICA; + else + return octave_syntax_detected ? LANG_OCTAVE : LANG_MATLAB; } #include diff --git a/test/expected_dir/mathematica1.m b/test/expected_dir/mathematica1.m index 43fcb2a..ea293c3 100644 --- a/test/expected_dir/mathematica1.m +++ b/test/expected_dir/mathematica1.m @@ -1466,9 +1466,9 @@ mathematica comment (* These are the arguments to createComment: mathematica comment - use vectorisation: False or True mathematica comment - use OpenCL: False or True mathematica comment - timelevels: 2 or 3 -mathematica comment (keep this at 3; this is better chosen with a run-time parameter) +mathematica comment ## (keep this at 3; this is better chosen with a run-time parameter) mathematica comment - matter: 0 or 1 -mathematica comment (matter seems cheap; it should be always enabled) +mathematica comment ## (matter seems cheap; it should be always enabled) mathematica comment - thorn base name mathematica comment *) mathematica blank diff --git a/test/expected_dir/octave1.m b/test/expected_dir/octave1.m index 311ca33..73977b0 100644 --- a/test/expected_dir/octave1.m +++ b/test/expected_dir/octave1.m @@ -4,7 +4,7 @@ octave comment ## 3-clause BSD license appended to this file. octave blank octave code function varargout = toledolu(LU) -octave comment ## -*- texinfo -*- +octave comment ## (*- texinfo -*) octave comment ## @deftypefn{Function File} {[@var{L}, @var{U}, @var{P}]} = toledolu(@var{A}) octave comment ## @deftypefnx{Function File} {[@var{L}, @var{U}]} = toledolu(@var{A}) octave comment ## @deftypefnx{Function File} {@var{LUP}} = toledolu(@var{A}) @@ -20,8 +20,8 @@ octave comment ## See the help for lu for details about the other calling forms. octave comment ## octave comment ## For the algorithm, see -octave comment ## @itemize -octave comment ## @item +octave comment ## (* @itemize *) +octave comment ## (* @item *) octave comment ## Toledo, Sivan. "Locality of reference in LU decomposition with octave comment ## partial pivoting," SIAM J. of Matrix Analysis and Applications, octave comment ## v18, n4, 1997. DOI: 10.1137/S0895479896297744 diff --git a/test/src_dir/mathematica1.m b/test/src_dir/mathematica1.m index 4bec755..6cccde6 100644 --- a/test/src_dir/mathematica1.m +++ b/test/src_dir/mathematica1.m @@ -1466,9 +1466,9 @@ which is not always the case (since we don't enforce it). On - use vectorisation: False or True - use OpenCL: False or True - timelevels: 2 or 3 - (keep this at 3; this is better chosen with a run-time parameter) + ## (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) + ## (matter seems cheap; it should be always enabled) - thorn base name *) diff --git a/test/src_dir/octave1.m b/test/src_dir/octave1.m index d72372d..4c795fb 100644 --- a/test/src_dir/octave1.m +++ b/test/src_dir/octave1.m @@ -4,7 +4,7 @@ ## 3-clause BSD license appended to this file. function varargout = toledolu(LU) - ## -*- texinfo -*- + ## (*- texinfo -*) ## @deftypefn{Function File} {[@var{L}, @var{U}, @var{P}]} = toledolu(@var{A}) ## @deftypefnx{Function File} {[@var{L}, @var{U}]} = toledolu(@var{A}) ## @deftypefnx{Function File} {@var{LUP}} = toledolu(@var{A}) @@ -20,8 +20,8 @@ ## See the help for lu for details about the other calling forms. ## ## For the algorithm, see - ## @itemize - ## @item + ## (* @itemize *) + ## (* @item *) ## Toledo, Sivan. "Locality of reference in LU decomposition with ## partial pivoting," SIAM J. of Matrix Analysis and Applications, ## v18, n4, 1997. DOI: 10.1137/S0895479896297744 From 86c728432709b84f5a65ad33822b4e661135255f Mon Sep 17 00:00:00 2001 From: Alex Johnson Date: Thu, 29 Jun 2017 16:32:14 +0530 Subject: [PATCH 3/3] OTWO-3061 Detect alternate mathematica extensions --- src/hash/extensions.gperf | 3 +++ test/unit/detector_test.h | 3 +++ 2 files changed, 6 insertions(+) diff --git a/src/hash/extensions.gperf b/src/hash/extensions.gperf index 0fb0e4a..fad70e7 100755 --- a/src/hash/extensions.gperf +++ b/src/hash/extensions.gperf @@ -133,6 +133,7 @@ mov, BINARY mp, LANG_METAPOST_WITH_TEX mp3, BINARY mpg, BINARY +mt, LANG_MATHEMATICA mustache, LANG_HTML mxml, LANG_MXML nix, LANG_NIX @@ -217,6 +218,8 @@ vhd, LANG_VHDL vhdl, LANG_VHDL vim, LANG_VIM wav, BINARY +wl, LANG_MATHEMATICA +wlt, LANG_MATHEMATICA xaml, LANG_XAML xls, BINARY xlw, BINARY diff --git a/test/unit/detector_test.h b/test/unit/detector_test.h index 64735a2..4719cbd 100755 --- a/test/unit/detector_test.h +++ b/test/unit/detector_test.h @@ -116,6 +116,9 @@ void test_detector_detect_polyglot() { ASSERT_DETECT(LANG_CPP, "uses_cpp_keywords.h"); ASSERT_DETECT(LANG_RUBY, "foo.rb"); ASSERT_DETECT(LANG_MAKE, "foo.mk"); + ASSERT_DETECT(LANG_MATHEMATICA, "foo.mt"); + ASSERT_DETECT(LANG_MATHEMATICA, "foo.wl"); + ASSERT_DETECT(LANG_MATHEMATICA, "foo.wlt"); ASSERT_DETECT(LANG_OBJECTIVE_C, "foo_objective_c.h"); ASSERT_DETECT(LANG_PHP, "upper_case_php"); ASSERT_DETECT(LANG_SMALLTALK, "example.st");