Correct initialisation of dtlapse

5 files changed, 75 insertions, 43 deletions

diff --git a/ML_BSSN/schedule.ccl b/ML_BSSN/schedule.ccl
index 6b870ab..f6fec8f 100644
--- a/ML_BSSN/schedule.ccl
+++ b/ML_BSSN/schedule.ccl
@@ -95,6 +95,8 @@ if (CCTK_EQUALS(my_initial_data, "ADMBase"))
   {
     LANG: C
     
+    SYNC: ML_dtlapse
+    SYNC: ML_dtshift
     SYNC: ML_Gamma
   } "ML_BSSN_convertFromADMBaseGamma"
 }
@@ -153,6 +155,12 @@ schedule ML_BSSN_convertToADMBase IN MoL_PostStep AFTER ML_BSSN_ApplyBCs AFTER M
 {
   LANG: C
   
+  SYNC: ADMBase::curv
+  SYNC: ADMBase::dtlapse
+  SYNC: ADMBase::dtshift
+  SYNC: ADMBase::lapse
+  SYNC: ADMBase::metric
+  SYNC: ADMBase::shift
 } "ML_BSSN_convertToADMBase"
 
 schedule ML_BSSN_constraints AT analysis
diff --git a/ML_BSSN/src/ML_BSSN_convertFromADMBase.c b/ML_BSSN/src/ML_BSSN_convertFromADMBase.c
index 599f013..7f156c5 100644
--- a/ML_BSSN/src/ML_BSSN_convertFromADMBase.c
+++ b/ML_BSSN/src/ML_BSSN_convertFromADMBase.c
@@ -101,19 +101,13 @@ void ML_BSSN_convertFromADMBase_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, C
     CCTK_REAL K11 = INITVALUE, K12 = INITVALUE, K13 = INITVALUE, K22 = INITVALUE, K23 = INITVALUE, K33 = INITVALUE;
     
     /* Declare local copies of grid functions */
-    CCTK_REAL AL = INITVALUE;
     CCTK_REAL alpL = INITVALUE;
     CCTK_REAL alphaL = INITVALUE;
     CCTK_REAL At11L = INITVALUE, At12L = INITVALUE, At13L = INITVALUE, At22L = INITVALUE, At23L = INITVALUE, At33L = INITVALUE;
-    CCTK_REAL B1L = INITVALUE, B2L = INITVALUE, B3L = INITVALUE;
     CCTK_REAL beta1L = INITVALUE, beta2L = INITVALUE, beta3L = INITVALUE;
     CCTK_REAL betaxL = INITVALUE;
     CCTK_REAL betayL = INITVALUE;
     CCTK_REAL betazL = INITVALUE;
-    CCTK_REAL dtalpL = INITVALUE;
-    CCTK_REAL dtbetaxL = INITVALUE;
-    CCTK_REAL dtbetayL = INITVALUE;
-    CCTK_REAL dtbetazL = INITVALUE;
     CCTK_REAL gt11L = INITVALUE, gt12L = INITVALUE, gt13L = INITVALUE, gt22L = INITVALUE, gt23L = INITVALUE, gt33L = INITVALUE;
     CCTK_REAL gxxL = INITVALUE;
     CCTK_REAL gxyL = INITVALUE;
@@ -138,10 +132,6 @@ void ML_BSSN_convertFromADMBase_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, C
     betaxL = betax[index];
     betayL = betay[index];
     betazL = betaz[index];
-    dtalpL = dtalp[index];
-    dtbetaxL = dtbetax[index];
-    dtbetayL = dtbetay[index];
-    dtbetazL = dtbetaz[index];
     gxxL = gxx[index];
     gxyL = gxy[index];
     gxzL = gxz[index];
@@ -236,23 +226,14 @@ void ML_BSSN_convertFromADMBase_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, C
     
     alphaL  =  alpL;
     
-    AL  =  dtalpL;
-    
     beta1L  =  betaxL;
     
     beta2L  =  betayL;
     
     beta3L  =  betazL;
     
-    B1L  =  dtbetaxL;
-    
-    B2L  =  dtbetayL;
-    
-    B3L  =  dtbetazL;
-    
     
     /* Copy local copies back to grid functions */
-    A[index] = AL;
     alpha[index] = alphaL;
     At11[index] = At11L;
     At12[index] = At12L;
@@ -260,9 +241,6 @@ void ML_BSSN_convertFromADMBase_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, C
     At22[index] = At22L;
     At23[index] = At23L;
     At33[index] = At33L;
-    B1[index] = B1L;
-    B2[index] = B2L;
-    B3[index] = B3L;
     beta1[index] = beta1L;
     beta2[index] = beta2L;
     beta3[index] = beta3L;
diff --git a/ML_BSSN/src/ML_BSSN_convertFromADMBaseGamma.c b/ML_BSSN/src/ML_BSSN_convertFromADMBaseGamma.c
index 3eb0c47..33930ba 100644
--- a/ML_BSSN/src/ML_BSSN_convertFromADMBaseGamma.c
+++ b/ML_BSSN/src/ML_BSSN_convertFromADMBaseGamma.c
@@ -101,11 +101,22 @@ void ML_BSSN_convertFromADMBaseGamma_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT fa
     CCTK_REAL gtu11 = INITVALUE, gtu21 = INITVALUE, gtu22 = INITVALUE, gtu31 = INITVALUE, gtu32 = INITVALUE, gtu33 = INITVALUE;
     
     /* Declare local copies of grid functions */
+    CCTK_REAL AL = INITVALUE;
+    CCTK_REAL alphaL = INITVALUE;
+    CCTK_REAL B1L = INITVALUE, B2L = INITVALUE, B3L = INITVALUE;
+    CCTK_REAL beta1L = INITVALUE, beta2L = INITVALUE, beta3L = INITVALUE;
+    CCTK_REAL dtalpL = INITVALUE;
+    CCTK_REAL dtbetaxL = INITVALUE;
+    CCTK_REAL dtbetayL = INITVALUE;
+    CCTK_REAL dtbetazL = INITVALUE;
     CCTK_REAL gt11L = INITVALUE, gt12L = INITVALUE, gt13L = INITVALUE, gt22L = INITVALUE, gt23L = INITVALUE, gt33L = INITVALUE;
     CCTK_REAL Xt1L = INITVALUE, Xt2L = INITVALUE, Xt3L = INITVALUE;
     /* Declare precomputed derivatives*/
     
     /* Declare derivatives */
+    CCTK_REAL PDstandardNth1alpha = INITVALUE;
+    CCTK_REAL PDstandardNth2alpha = INITVALUE;
+    CCTK_REAL PDstandardNth3alpha = INITVALUE;
     CCTK_REAL PDstandardNth1gt11 = INITVALUE;
     CCTK_REAL PDstandardNth2gt11 = INITVALUE;
     CCTK_REAL PDstandardNth3gt11 = INITVALUE;
@@ -126,6 +137,14 @@ void ML_BSSN_convertFromADMBaseGamma_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT fa
     CCTK_REAL PDstandardNth3gt33 = INITVALUE;
     
     /* Assign local copies of grid functions */
+    alphaL = alpha[index];
+    beta1L = beta1[index];
+    beta2L = beta2[index];
+    beta3L = beta3[index];
+    dtalpL = dtalp[index];
+    dtbetaxL = dtbetax[index];
+    dtbetayL = dtbetay[index];
+    dtbetazL = dtbetaz[index];
     gt11L = gt11[index];
     gt12L = gt12[index];
     gt13L = gt13[index];
@@ -138,6 +157,9 @@ void ML_BSSN_convertFromADMBaseGamma_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT fa
     /* Include user supplied include files */
     
     /* Precompute derivatives (new style) */
+    PDstandardNth1alpha = PDstandardNth1(alpha, i, j, k);
+    PDstandardNth2alpha = PDstandardNth2(alpha, i, j, k);
+    PDstandardNth3alpha = PDstandardNth3(alpha, i, j, k);
     PDstandardNth1gt11 = PDstandardNth1(gt11, i, j, k);
     PDstandardNth2gt11 = PDstandardNth2(gt11, i, j, k);
     PDstandardNth3gt11 = PDstandardNth3(gt11, i, j, k);
@@ -234,8 +256,21 @@ void ML_BSSN_convertFromADMBaseGamma_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT fa
     
     Xt3L  =  Gt311*gtu11 + Gt322*gtu22 + 2*(Gt312*gtu21 + Gt313*gtu31 + Gt323*gtu32) + Gt333*gtu33;
     
+    AL  =  (-dtalpL + beta1L*PDstandardNth1alpha + beta2L*PDstandardNth2alpha + beta3L*PDstandardNth3alpha)*INV(harmonicF)*
+        pow(alphaL,-harmonicN);
+    
+    B1L  =  dtbetaxL*ShiftGammaCoeff*INV(1.e-100 + SQR(ShiftGammaCoeff))*pow(alphaL,-ShiftAlphaPower);
+    
+    B2L  =  dtbetayL*ShiftGammaCoeff*INV(1.e-100 + SQR(ShiftGammaCoeff))*pow(alphaL,-ShiftAlphaPower);
+    
+    B3L  =  dtbetazL*ShiftGammaCoeff*INV(1.e-100 + SQR(ShiftGammaCoeff))*pow(alphaL,-ShiftAlphaPower);
+    
     
     /* Copy local copies back to grid functions */
+    A[index] = AL;
+    B1[index] = B1L;
+    B2[index] = B2L;
+    B3[index] = B3L;
     Xt1[index] = Xt1L;
     Xt2[index] = Xt2L;
     Xt3[index] = Xt3L;
diff --git a/ML_BSSN/src/ML_BSSN_convertToADMBase.c b/ML_BSSN/src/ML_BSSN_convertToADMBase.c
index f51b779..0586578 100644
--- a/ML_BSSN/src/ML_BSSN_convertToADMBase.c
+++ b/ML_BSSN/src/ML_BSSN_convertToADMBase.c
@@ -130,6 +130,9 @@ void ML_BSSN_convertToADMBase_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCT
     /* Declare precomputed derivatives*/
     
     /* Declare derivatives */
+    CCTK_REAL PDstandardNth1alpha = INITVALUE;
+    CCTK_REAL PDstandardNth2alpha = INITVALUE;
+    CCTK_REAL PDstandardNth3alpha = INITVALUE;
     
     /* Assign local copies of grid functions */
     AL = A[index];
@@ -160,6 +163,9 @@ void ML_BSSN_convertToADMBase_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCT
     /* Include user supplied include files */
     
     /* Precompute derivatives (new style) */
+    PDstandardNth1alpha = PDstandardNth1(alpha, i, j, k);
+    PDstandardNth2alpha = PDstandardNth2(alpha, i, j, k);
+    PDstandardNth3alpha = PDstandardNth3(alpha, i, j, k);
     
     /* Precompute derivatives (old style) */
     
@@ -216,19 +222,20 @@ void ML_BSSN_convertToADMBase_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCT
     
     alpL  =  alphaL;
     
-    dtalpL  =  AL;
-    
     betaxL  =  beta1L;
     
     betayL  =  beta2L;
     
     betazL  =  beta3L;
     
-    dtbetaxL  =  B1L;
+    dtalpL  =  beta1L*PDstandardNth1alpha + beta2L*PDstandardNth2alpha + beta3L*PDstandardNth3alpha - 
+        AL*harmonicF*pow(alphaL,harmonicN);
+    
+    dtbetaxL  =  B1L*ShiftGammaCoeff*pow(alphaL,ShiftAlphaPower);
     
-    dtbetayL  =  B2L;
+    dtbetayL  =  B2L*ShiftGammaCoeff*pow(alphaL,ShiftAlphaPower);
     
-    dtbetazL  =  B3L;
+    dtbetazL  =  B3L*ShiftGammaCoeff*pow(alphaL,ShiftAlphaPower);
     
     
     /* Copy local copies back to grid functions */
@@ -263,5 +270,5 @@ void ML_BSSN_convertToADMBase(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_PARAMETERS
   
-  GenericFD_LoopOverEverything(cctkGH, &ML_BSSN_convertToADMBase_Body);
+  GenericFD_LoopOverInterior(cctkGH, &ML_BSSN_convertToADMBase_Body);
 }
diff --git a/m/McLachlan.m b/m/McLachlan.m
index 6973069..d27041d 100644
--- a/m/McLachlan.m
+++ b/m/McLachlan.m
@@ -258,16 +258,10 @@ convertFromADMBaseCalcBSSN =
     At[la,lb] -> em4phi (K[la,lb] - (1/3) g[la,lb] trK),
     
     alpha -> alp,
-    (* TODO: this is wrong *)
-    A     -> dtalp,
     
     beta1 -> betax,
     beta2 -> betay,
-    beta3 -> betaz,
-    (* TODO: this is wrong *)
-    B1    -> dtbetax,
-    B2    -> dtbetay,
-    B3    -> dtbetaz
+    beta3 -> betaz
   }
 }
 
@@ -284,7 +278,18 @@ convertFromADMBaseCalcBSSNGamma =
     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]
+    Xt[ua] -> gtu[ub,uc] Gt[ua,lb,lc],
+    
+    (* TODO: check this *)
+    A -> - (dtalp - Lie[alpha, beta]) / (harmonicF alpha^harmonicN),
+    
+    (* TODO: check this *)
+    (* B1 -> dtbetax / (ShiftGammaCoeff alpha^ShiftAlphaPower) *)
+    (* B2 -> dtbetay / (ShiftGammaCoeff alpha^ShiftAlphaPower) *)
+    (* B3 -> dtbetaz / (ShiftGammaCoeff alpha^ShiftAlphaPower) *)
+    B1 -> ShiftGammaCoeff dtbetax / ((ShiftGammaCoeff^2 + 1.0*^-100) alpha^ShiftAlphaPower),
+    B2 -> ShiftGammaCoeff dtbetay / ((ShiftGammaCoeff^2 + 1.0*^-100) alpha^ShiftAlphaPower),
+    B3 -> ShiftGammaCoeff dtbetaz / ((ShiftGammaCoeff^2 + 1.0*^-100) alpha^ShiftAlphaPower)
   }
 }
 
@@ -326,6 +331,7 @@ convertToADMBaseCalcBSSN =
 {
   Name -> "ML_BSSN_convertToADMBase",
   Schedule -> {"IN MoL_PostStep AFTER ML_BSSN_ApplyBCs AFTER ML_BSSN_enforce"},
+  Where -> Interior,
   Shorthands -> {e4phi, g[la,lb], K[la,lb]},
   Equations -> 
   {
@@ -345,16 +351,14 @@ convertToADMBaseCalcBSSN =
     kyz      -> K23,
     kzz      -> K33,
     alp      -> alpha,
-    (* TODO: this is wrong *)
-    (* TODO: rename dtalp->A, dtbeta->B *)
-    dtalp    -> A,
     betax    -> beta1,
     betay    -> beta2,
     betaz    -> beta3,
-    (* TODO: this is wrong *)
-    dtbetax  -> B1,
-    dtbetay  -> B2,
-    dtbetaz  -> B3
+    (* see RHS *)
+    dtalp    -> - harmonicF alpha^harmonicN A + Lie[alpha, beta],
+    dtbetax  -> ShiftGammaCoeff alpha^ShiftAlphaPower B1,
+    dtbetay  -> ShiftGammaCoeff alpha^ShiftAlphaPower B2,
+    dtbetaz  -> ShiftGammaCoeff alpha^ShiftAlphaPower B3
   }
 }