Add gauge conditions

4 files changed, 104 insertions, 52 deletions

diff --git a/ML_BSSN/param.ccl b/ML_BSSN/param.ccl
index 918049a..1430b55 100644
--- a/ML_BSSN/param.ccl
+++ b/ML_BSSN/param.ccl
@@ -29,6 +29,42 @@ CCTK_REAL verbose "verbose"
   "*:*" :: "no restrictions"
 } 0
 
+restricted:
+CCTK_REAL harmonicF "harmonicF"
+{
+  "*:*" :: "no restrictions"
+} 1
+
+restricted:
+CCTK_REAL AlphaDriver "AlphaDriver"
+{
+  "*:*" :: "no restrictions"
+} 0
+
+restricted:
+CCTK_REAL ShiftGammaCoeff "ShiftGammaCoeff"
+{
+  "*:*" :: "no restrictions"
+} 0
+
+restricted:
+CCTK_REAL BetaDriver "BetaDriver"
+{
+  "*:*" :: "no restrictions"
+} 0
+
+restricted:
+CCTK_INT harmonicN "harmonicN"
+{
+  *:* :: "no restrictions"
+} 1
+
+restricted:
+CCTK_INT ShiftAlphaPower "ShiftAlphaPower"
+{
+  *:* :: "no restrictions"
+} 0
+
 private:
 KEYWORD my_initial_data "my_initial_data"
 {
diff --git a/ML_BSSN/src/ML_BSSN_RHS.c b/ML_BSSN/src/ML_BSSN_RHS.c
index 930a860..a11173b 100644
--- a/ML_BSSN/src/ML_BSSN_RHS.c
+++ b/ML_BSSN/src/ML_BSSN_RHS.c
@@ -212,18 +212,6 @@ void ML_BSSN_RHS_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL normal
     CCTK_REAL PDstandard4th23beta3 = INITVALUE;
     CCTK_REAL PDstandard4th31beta3 = INITVALUE;
     CCTK_REAL PDstandard4th32beta3 = INITVALUE;
-    CCTK_REAL PDstandard4th1dtalpha = INITVALUE;
-    CCTK_REAL PDstandard4th2dtalpha = INITVALUE;
-    CCTK_REAL PDstandard4th3dtalpha = INITVALUE;
-    CCTK_REAL PDstandard4th1dtbeta1 = INITVALUE;
-    CCTK_REAL PDstandard4th2dtbeta1 = INITVALUE;
-    CCTK_REAL PDstandard4th3dtbeta1 = INITVALUE;
-    CCTK_REAL PDstandard4th1dtbeta2 = INITVALUE;
-    CCTK_REAL PDstandard4th2dtbeta2 = INITVALUE;
-    CCTK_REAL PDstandard4th3dtbeta2 = INITVALUE;
-    CCTK_REAL PDstandard4th1dtbeta3 = INITVALUE;
-    CCTK_REAL PDstandard4th2dtbeta3 = INITVALUE;
-    CCTK_REAL PDstandard4th3dtbeta3 = INITVALUE;
     CCTK_REAL PDstandard4th1gt11 = INITVALUE;
     CCTK_REAL PDstandard4th2gt11 = INITVALUE;
     CCTK_REAL PDstandard4th3gt11 = INITVALUE;
@@ -344,9 +332,13 @@ void ML_BSSN_RHS_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL normal
     gt33L = gt33[index];
     phiL = phi[index];
     trKL = trK[index];
+    trKrhsL = trKrhs[index];
     Xt1L = Xt1[index];
+    Xt1rhsL = Xt1rhs[index];
     Xt2L = Xt2[index];
+    Xt2rhsL = Xt2rhs[index];
     Xt3L = Xt3[index];
+    Xt3rhsL = Xt3rhs[index];
     
     /* Assign local copies of subblock grid functions */
     
@@ -407,18 +399,6 @@ void ML_BSSN_RHS_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL normal
     PDstandard4th12beta3 = PDstandard4th12(beta3, i, j, k);
     PDstandard4th13beta3 = PDstandard4th13(beta3, i, j, k);
     PDstandard4th23beta3 = PDstandard4th23(beta3, i, j, k);
-    PDstandard4th1dtalpha = PDstandard4th1(dtalpha, i, j, k);
-    PDstandard4th2dtalpha = PDstandard4th2(dtalpha, i, j, k);
-    PDstandard4th3dtalpha = PDstandard4th3(dtalpha, i, j, k);
-    PDstandard4th1dtbeta1 = PDstandard4th1(dtbeta1, i, j, k);
-    PDstandard4th2dtbeta1 = PDstandard4th2(dtbeta1, i, j, k);
-    PDstandard4th3dtbeta1 = PDstandard4th3(dtbeta1, i, j, k);
-    PDstandard4th1dtbeta2 = PDstandard4th1(dtbeta2, i, j, k);
-    PDstandard4th2dtbeta2 = PDstandard4th2(dtbeta2, i, j, k);
-    PDstandard4th3dtbeta2 = PDstandard4th3(dtbeta2, i, j, k);
-    PDstandard4th1dtbeta3 = PDstandard4th1(dtbeta3, i, j, k);
-    PDstandard4th2dtbeta3 = PDstandard4th2(dtbeta3, i, j, k);
-    PDstandard4th3dtbeta3 = PDstandard4th3(dtbeta3, i, j, k);
     PDstandard4th1gt11 = PDstandard4th1(gt11, i, j, k);
     PDstandard4th2gt11 = PDstandard4th2(gt11, i, j, k);
     PDstandard4th3gt11 = PDstandard4th3(gt11, i, j, k);
@@ -1448,24 +1428,22 @@ void ML_BSSN_RHS_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL normal
            2*alphaL*g33*gu23*R32 - 3*alphaL*R33 + alphaL*g33*gu33*R33) + 
         alphaL*(-2*(At31L*Atm13 + At32L*Atm23 + At33L*Atm33) + At33L*trKL);
     
-    alpharhsL  =  beta1L*PDstandard4th1alpha + beta2L*PDstandard4th2alpha + beta3L*PDstandard4th3alpha;
+    alpharhsL  =  beta1L*PDstandard4th1alpha + beta2L*PDstandard4th2alpha + beta3L*PDstandard4th3alpha - 
+        dtalphaL*harmonicF*pow(alphaL,harmonicN);
     
-    dtalpharhsL  =  beta1L*PDstandard4th1dtalpha + beta2L*PDstandard4th2dtalpha + beta3L*PDstandard4th3dtalpha;
+    dtalpharhsL  =  -(AlphaDriver*dtalphaL) + trKrhsL;
     
-    beta1rhsL  =  0;
+    beta1rhsL  =  dtbeta1L*ShiftGammaCoeff*pow(alphaL,ShiftAlphaPower);
     
-    beta2rhsL  =  0;
+    beta2rhsL  =  dtbeta2L*ShiftGammaCoeff*pow(alphaL,ShiftAlphaPower);
     
-    beta3rhsL  =  0;
+    beta3rhsL  =  dtbeta3L*ShiftGammaCoeff*pow(alphaL,ShiftAlphaPower);
     
-    dtbeta1rhsL  =  -(dtbeta1L*PDstandard4th1beta1) + beta1L*PDstandard4th1dtbeta1 - dtbeta2L*PDstandard4th2beta1 + 
-        beta2L*PDstandard4th2dtbeta1 - dtbeta3L*PDstandard4th3beta1 + beta3L*PDstandard4th3dtbeta1;
+    dtbeta1rhsL  =  -(BetaDriver*dtbeta1L) + Xt1rhsL;
     
-    dtbeta2rhsL  =  -(dtbeta1L*PDstandard4th1beta2) + beta1L*PDstandard4th1dtbeta2 - dtbeta2L*PDstandard4th2beta2 + 
-        beta2L*PDstandard4th2dtbeta2 - dtbeta3L*PDstandard4th3beta2 + beta3L*PDstandard4th3dtbeta2;
+    dtbeta2rhsL  =  -(BetaDriver*dtbeta2L) + Xt2rhsL;
     
-    dtbeta3rhsL  =  -(dtbeta1L*PDstandard4th1beta3) + beta1L*PDstandard4th1dtbeta3 - dtbeta2L*PDstandard4th2beta3 + 
-        beta2L*PDstandard4th2dtbeta3 - dtbeta3L*PDstandard4th3beta3 + beta3L*PDstandard4th3dtbeta3;
+    dtbeta3rhsL  =  -(BetaDriver*dtbeta3L) + Xt3rhsL;
     
     
     /* Copy local copies back to grid functions */
diff --git a/m/McLachlan.m b/m/McLachlan.m
index 559da75..cc2889a 100644
--- a/m/McLachlan.m
+++ b/m/McLachlan.m
@@ -462,14 +462,13 @@ evolCalcBSSN =
                       + alpha (trK At[la,lb] - 2 At[la,lc] Atm[uc,lb])
                       + Lie[At[la,lb], beta] - (2/3) At[la,lb] PD[beta[uc],lc],
     
-    dot[alpha]      -> (* TODO *)
+    (* dot[alpha] -> - harmonicF alpha^harmonicN trK, *)
+    dot[alpha]      -> - harmonicF alpha^harmonicN dtalpha
                        + Lie[alpha, beta],
-    dot[dtalpha]    -> (* TODO *)
-                       + Lie[dtalpha, beta],
-    dot[beta[ua]]   -> (* TODO *)
-                       + Lie[beta[ua], beta],
-    dot[dtbeta[ua]] -> (* TODO *)
-                       + Lie[dtbeta[ua], beta]
+    dot[dtalpha]    -> dot[trK] - AlphaDriver dtalpha,
+    (* dot[beta[ua]] -> eta Xt[ua], *)
+    dot[beta[ua]]   -> ShiftGammaCoeff alpha^ShiftAlphaPower dtbeta[ua],
+    dot[dtbeta[ua]] -> dot[Xt[ua]] - BetaDriver dtbeta[ua]
   }
 }
 
@@ -595,15 +594,49 @@ inheritedImplementations = {"ADMBase"};
 inheritedKeywordParameters = {"ADMBase::initial_data"};
 
 keywordParameters =
-{{
-  Name -> "my_initial_data",
-  (* Visibility -> "restricted", *)
-  (* Description -> "ddd", *)
-  AllowedValues -> {"ADMBase", "Minkowski"},
-  Default -> "ADMBase"
-}};
+{
+  {
+    Name -> "my_initial_data",
+    (* Visibility -> "restricted", *)
+    (* Description -> "ddd", *)
+    AllowedValues -> {"ADMBase", "Minkowski"},
+    Default -> "ADMBase"
+  }
+};
 
-realParameters = {};
+intParameters =
+{
+  {
+    Name -> harmonicN,
+    Description -> "d/dt alpha = - f alpha^n K  (harmonic=1, 1+log=0)",
+    Default -> 1
+  },
+  {
+    Name -> ShiftAlphaPower,
+    Default -> 0
+  }
+};
+
+realParameters =
+{
+  {
+    Name -> harmonicF,
+    Description -> "d/dt alpha = - f alpha^n K   (harmonic=1, 1+log=2)",
+    Default -> 1
+  },
+  {
+    Name -> AlphaDriver,
+    Default -> 0
+  },
+  {
+    Name -> ShiftGammaCoeff,
+    Default -> 0
+  },
+  {
+    Name -> BetaDriver,
+    Default -> 0
+  }
+};
 
 (******************************************************************************)
 (* Construct the thorns *)
@@ -624,8 +657,7 @@ CreateKrancThornTT [groups, ".", "ML_ADM",
   PartialDerivatives -> derivatives,
   InheritedImplementations -> inheritedImplementations,
   InheritedKeywordParameters -> inheritedKeywordParameters,
-  KeywordParameters -> keywordParameters,
-  RealParameters -> realParameters
+  KeywordParameters -> keywordParameters
 ];
 
 calculationsBSSN = 
@@ -646,5 +678,6 @@ CreateKrancThornTT [groupsBSSN, ".", "ML_BSSN",
   InheritedImplementations -> inheritedImplementations,
   InheritedKeywordParameters -> inheritedKeywordParameters,
   KeywordParameters -> keywordParameters,
+  IntParameters -> intParameters,
   RealParameters -> realParameters
 ];
diff --git a/par/m_bssn.par b/par/m_bssn.par
index f55771e..5533f4d 100644
--- a/par/m_bssn.par
+++ b/par/m_bssn.par
@@ -29,6 +29,11 @@ ADMBase::initial_dtshift = "zero"
 
 ML_BSSN::my_initial_data = "Minkowski"
 
+ML_BSSN::harmonic_n      = 0      # 1+log
+ML_BSSN::harmonic_f      = 2.0    # 1+log
+ML_BSSN::ShiftGammaCoeff = 0.75
+ML_BSSN::BetaDriver      = 0.5
+
 
 
 IO::out_dir = $parfile