From 71894c54a1d778df2b27277c004157aa6be93da7 Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Mon, 2 Jan 2012 11:58:25 -0500 Subject: Don't replace Pi by its numeric equivalent -- Kranc does this already --- m/McLachlanW.m | 12 +++++------- m/McLachlan_ADM.m | 2 -- m/McLachlan_ADMConstraints.m | 8 ++------ m/McLachlan_ADMQuantities.m | 10 ++++------ m/McLachlan_BSSN.m | 12 +++++------- m/McLachlan_WeylScalars.m | 6 ------ 6 files changed, 16 insertions(+), 34 deletions(-) (limited to 'm') diff --git a/m/McLachlanW.m b/m/McLachlanW.m index 4883dde..e8a329b 100644 --- a/m/McLachlanW.m +++ b/m/McLachlanW.m @@ -107,8 +107,6 @@ ddetgtExpr[la_] = Sum [D[Det[MatrixOfComponents[gt[la, lb]]], X] PD[X, la], {X, Union[Flatten[MatrixOfComponents[gt[la, lb]]]]}]; -pi = N[Pi,40]; - (******************************************************************************) (* Groups *) (******************************************************************************) @@ -481,14 +479,14 @@ evolCalcBSSNW = - Xtn[uj] PD[beta[ui],lj] + (2/3) Xtn[ui] PD[beta[uj],lj] (* Equation (4.28) in Baumgarte & Shapiro (Phys. Rept. 376 (2003) 41-131) *) - + addMatter (- 16 pi alpha gtu[ui,uj] S[lj]), + + addMatter (- 16 Pi alpha gtu[ui,uj] S[lj]), (* PRD 62, 044034 (2000), eqn. (11) *) dot[trK] -> - gu[ua,ub] CD[alpha,la,lb] + alpha (Atm[ua,lb] Atm[ub,la] + (1/3) trK^2) + ( betam[ua] PDm[trK,la] + betap[ua] PDp[trK,la] ) (* Equation (4.21) in Baumgarte & Shapiro (Phys. Rept. 376 (2003) 41-131) *) - + addMatter (4 pi alpha (rho + trS)), + + addMatter (4 Pi alpha (rho + trS)), (* PRD 62, 044034 (2000), eqn. (12) *) (* TODO: use Hamiltonian constraint to make tracefree *) @@ -501,7 +499,7 @@ evolCalcBSSNW = + 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 (- W2 alpha 8 pi + + addMatter (- W2 alpha 8 Pi (T[la,lb] - (1/3) g[la,lb] trS)), (* dot[alpha] -> - harmonicF alpha^harmonicN trK, *) @@ -682,7 +680,7 @@ constraintsCalcBSSNW = (* 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, + 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] @@ -692,7 +690,7 @@ constraintsCalcBSSNW = M[li] -> + gtu[uj,uk] (CDt[At[li,lj],lk] + 6 At[li,lj] pdphi[lk]) - (2/3) PD[trK,li] - - addMatter 8 pi S[li], + - addMatter 8 Pi S[li], (* TODO: use PRD 67, 084023 (2003), eqn. (20) *) (* det gamma-tilde *) diff --git a/m/McLachlan_ADM.m b/m/McLachlan_ADM.m index 92f0590..94ab2ae 100644 --- a/m/McLachlan_ADM.m +++ b/m/McLachlan_ADM.m @@ -82,8 +82,6 @@ ddetgExpr[la_] = Sum [D[Det[MatrixOfComponents[g[la, lb]]], X] PD[X, la], {X, Union[Flatten[MatrixOfComponents[g[la, lb]]]]}]; -pi = N[Pi,40]; - (******************************************************************************) (* Groups *) (******************************************************************************) diff --git a/m/McLachlan_ADMConstraints.m b/m/McLachlan_ADMConstraints.m index 19346e7..eaddae6 100644 --- a/m/McLachlan_ADMConstraints.m +++ b/m/McLachlan_ADMConstraints.m @@ -82,14 +82,10 @@ ddetgExpr[la_] = Sum [D[Det[MatrixOfComponents[g[la, lb]]], X] PD[X, la], {X, Union[Flatten[MatrixOfComponents[g[la, lb]]]]}]; -pi = N[Pi,40]; - (******************************************************************************) (* Groups *) (******************************************************************************) -SetGroupTimelevels[g_,tl_] = Join[g, {Timelevels -> tl}]; - evolvedGroups = {}; evaluatedGroups = {SetGroupName [CreateGroupFromTensor [H ], prefix <> "Ham"], @@ -163,9 +159,9 @@ ADMConstraintsCalc = (* ADM constraints *) H -> + trR - Km[ua,lb] Km[ub,la] + trK^2 - - addMatter 16 pi rho, + - addMatter 16 Pi rho, M[la] -> + gu[ub,uc] (CD[K[lc,la], lb] - CD[K[lc,lb], la]) - - addMatter 8 pi S[la] + - addMatter 8 Pi S[la] } }; diff --git a/m/McLachlan_ADMQuantities.m b/m/McLachlan_ADMQuantities.m index 003a9f8..2f78d56 100644 --- a/m/McLachlan_ADMQuantities.m +++ b/m/McLachlan_ADMQuantities.m @@ -149,8 +149,6 @@ T11=eTxx; T12=eTxy; T22=eTyy; T13=eTxz; T23=eTyz; T33=eTzz; detgtExpr = Det [MatrixOfComponents [gt[la,lb]]]; -pi = N[Pi,40]; - (******************************************************************************) (* Groups *) (******************************************************************************) @@ -257,16 +255,16 @@ ADMQuantitiesCalc = (* ADM quantities *) (* See PRD 66, 084026 (2002) *) - Madm -> 1/(16 pi) - (+ ephi^5 (+ 16 pi addMatter rho + Madm -> 1/(16 Pi) + (+ ephi^5 (+ 16 Pi addMatter rho + Atm[ua,lb] Atm[ub,la] - 2/3 trK^2) - gtu[ua,ub] Gt[uc,la,ld] Gtlu[lb,lc,ud] + (1 - ephi) trRt), - Jadm[li] -> 1/(16 pi) Eps[li,lj,uk] ephi^6 + Jadm[li] -> 1/(16 Pi) Eps[li,lj,uk] ephi^6 (+ 2 Atm[uj,lk] - + 16 pi x[uj] S[lk] + + 16 Pi x[uj] S[lk] + 4/3 x[uj] PD[trK,lk] - x[uj] dgtu[ul,um,lk] At[ll,lm]) } diff --git a/m/McLachlan_BSSN.m b/m/McLachlan_BSSN.m index 107ecd4..802ae57 100644 --- a/m/McLachlan_BSSN.m +++ b/m/McLachlan_BSSN.m @@ -159,8 +159,6 @@ T11=eTxx; T12=eTxy; T22=eTyy; T13=eTxz; T23=eTyz; T33=eTzz; (* Expressions *) (******************************************************************************) -pi = N[Pi,40]; - detgExpr = Det [MatrixOfComponents [g [la,lb]]]; ddetgExpr[la_] = Sum [D[Det[MatrixOfComponents[g[la, lb]]], X] PD[X, la], @@ -613,7 +611,7 @@ evolCalc = - Xtn[uj] PD[beta[ui],lj] + (2/3) Xtn[ui] PD[beta[uj],lj] (* Equation (4.28) in Baumgarte & Shapiro (Phys. Rept. 376 (2003) 41-131) *) - + addMatter (- 16 pi alpha gtu[ui,uj] S[lj]), + + addMatter (- 16 Pi alpha gtu[ui,uj] S[lj]), dot[Xt[ui]] -> dotXt[ui], (* PRD 62, 044034 (2000), eqn. (11) *) @@ -622,7 +620,7 @@ evolCalc = - Xtn[ua] PD[alpha,la] ) + alpha (Atm[ua,lb] Atm[ub,la] + (1/3) trK^2) (* Equation (4.21) in Baumgarte & Shapiro (Phys. Rept. 376 (2003) 41-131) *) - + addMatter (4 pi alpha (rho + trS)), + + addMatter (4 Pi alpha (rho + trS)), dot[trK] -> dottrK, (* PRD 62, 044034 (2000), eqn. (12) *) @@ -636,7 +634,7 @@ evolCalc = + 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 + + addMatter (- em4phi alpha 8 Pi (T[la,lb] - (1/3) g[la,lb] trS)), (* dot[alpha] -> - harmonicF alpha^harmonicN trK, *) @@ -1018,7 +1016,7 @@ constraintsCalc = (* 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, + 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] @@ -1028,7 +1026,7 @@ constraintsCalc = 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], + - addMatter 8 Pi S[li], (* TODO: use PRD 67, 084023 (2003), eqn. (20) *) (* det gamma-tilde *) diff --git a/m/McLachlan_WeylScalars.m b/m/McLachlan_WeylScalars.m index 5e8d805..df10897 100644 --- a/m/McLachlan_WeylScalars.m +++ b/m/McLachlan_WeylScalars.m @@ -104,12 +104,6 @@ Map [DefineTensor, betax, betay, betaz, dtbetax, dtbetay, dtbetaz}]; -(******************************************************************************) -(* Expressions *) -(******************************************************************************) - -pi = N[Pi,40]; - (******************************************************************************) (* Groups *) (******************************************************************************) -- cgit v1.2.3 From cea0b49075b30965932112399b8531292c5f09be Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Mon, 2 Jan 2012 11:58:59 -0500 Subject: Remove unused function definition SetGroupTimelevels --- m/McLachlan_BSSN.m | 2 -- 1 file changed, 2 deletions(-) (limited to 'm') diff --git a/m/McLachlan_BSSN.m b/m/McLachlan_BSSN.m index 802ae57..7ff370c 100644 --- a/m/McLachlan_BSSN.m +++ b/m/McLachlan_BSSN.m @@ -178,8 +178,6 @@ thetaExpr = Min [Exp [1 - r / SpatialShiftGammaCoeffRadius], 1]; (* Groups *) (******************************************************************************) -SetGroupTimelevels[g_,tl_] = Join[g, {Timelevels -> tl}]; - evolvedGroups = {SetGroupName [CreateGroupFromTensor [phi ], prefix <> "log_confac"], SetGroupName [CreateGroupFromTensor [gt[la,lb]], prefix <> "metric" ], -- cgit v1.2.3 From 6cd8c6e35730454e4091ca1fc11c5ded823d6448 Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Mon, 2 Jan 2012 12:02:53 -0500 Subject: Correct definition of dot[A] advection terms There terms were implemented wrong, probably when the advection terms were split out of the main RHS routine. Without this correction, Kerr-Schild is not manifestly stationary any more. Some of these terms are in code paths that are (to my knowledge) unused, and so it is not clear what to do. I've tried to make these code paths as similar to the others as possible. --- m/McLachlan_BSSN.m | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) (limited to 'm') diff --git a/m/McLachlan_BSSN.m b/m/McLachlan_BSSN.m index 7ff370c..4df6573 100644 --- a/m/McLachlan_BSSN.m +++ b/m/McLachlan_BSSN.m @@ -648,10 +648,9 @@ evolCalc = *) dot[alpha] -> - harmonicF alpha^harmonicN (+ LapseACoeff A - + (1 - LapseACoeff) trK), + + (1 - LapseACoeff) (trK + AlphaDriver (alpha - 1))), dot[A] -> + LapseACoeff (dottrK - AlphaDriver A), - eta -> etaExpr, theta -> thetaExpr, @@ -702,16 +701,21 @@ advectCalc = 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[alpha] -> dot[alpha] + + LapseAdvectionCoeff Upwind[beta[ua], alpha, la], - dot[A] -> dot[A] + LapseAdvectionCoeff Upwind[beta[ua], A, la], + dot[A] -> dot[A] + + LapseAdvectionCoeff Upwind[beta[ua], A, la] + + (LapseACoeff - LapseAdvectionCoeff) + Upwind[beta[ua], trK, la], - dot[beta[ua]] -> dot[beta[ua]] + ShiftAdvectionCoeff Upwind[beta[ub], beta[ua], lb], + dot[beta[ua]] -> dot[beta[ua]] + + ShiftAdvectionCoeff Upwind[beta[ub], beta[ua], lb], dot[B[ua]] -> dot[B[ua]] - + ShiftBCoeff Upwind[beta[uj], Xt[ua], lj] (* take care *) + ShiftAdvectionCoeff Upwind[beta[ub], B[ua], lb] - - ShiftAdvectionCoeff Upwind[beta[ub], Xt[ua], lb] + + (ShiftBCoeff - ShiftAdvectionCoeff) + Upwind[beta[ub], Xt[ua], lb] } }; -- cgit v1.2.3 From 8cb581f64d8116593bbc0b046571ac17e932b612 Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Wed, 4 Jan 2012 12:20:56 -0500 Subject: Apply advection terms to A or B^i only if they are evolved --- m/McLachlan_BSSN.m | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) (limited to 'm') diff --git a/m/McLachlan_BSSN.m b/m/McLachlan_BSSN.m index 4df6573..a8f6984 100644 --- a/m/McLachlan_BSSN.m +++ b/m/McLachlan_BSSN.m @@ -705,17 +705,18 @@ advectCalc = + LapseAdvectionCoeff Upwind[beta[ua], alpha, la], dot[A] -> dot[A] - + LapseAdvectionCoeff Upwind[beta[ua], A, la] - + (LapseACoeff - LapseAdvectionCoeff) - Upwind[beta[ua], trK, la], + + 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]] - + ShiftAdvectionCoeff Upwind[beta[ub], B[ua], lb] - + (ShiftBCoeff - ShiftAdvectionCoeff) - Upwind[beta[ub], Xt[ua], lb] + + ShiftBCoeff ( + + ShiftAdvectionCoeff Upwind[beta[ub], B[ua], lb] + + ((1 - ShiftAdvectionCoeff) + Upwind[beta[ub], Xt[ua], lb])) } }; -- cgit v1.2.3 From 45ec267196a8f999bd0bc0b72a18ce2bbcdbf4af Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Sun, 5 Feb 2012 21:54:48 -0500 Subject: Correct etaExpr to avoid (inconsequential) division by zero --- m/McLachlan_BSSN.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'm') diff --git a/m/McLachlan_BSSN.m b/m/McLachlan_BSSN.m index a8f6984..2d4eeb3 100644 --- a/m/McLachlan_BSSN.m +++ b/m/McLachlan_BSSN.m @@ -169,7 +169,7 @@ ddetgtExpr[la_] = Sum [D[Det[MatrixOfComponents[gt[la, lb]]], X] PD[X, la], {X, Union[Flatten[MatrixOfComponents[gt[la, lb]]]]}]; -etaExpr = Min [SpatialBetaDriverRadius / r, 1]; +etaExpr = SpatialBetaDriverRadius / Max [r, SpatialBetaDriverRadius]; thetaExpr = Min [Exp [1 - r / SpatialShiftGammaCoeffRadius], 1]; -- cgit v1.2.3