Add SpatialShiftGammaCoeffRadius to reduce ShiftGammaCoeff outside of a certain radius

6 files changed, 72 insertions, 38 deletions

diff --git a/ML_BSSN/param.ccl b/ML_BSSN/param.ccl
index 8fb0583..188e0f1 100644
--- a/ML_BSSN/param.ccl
+++ b/ML_BSSN/param.ccl
@@ -103,6 +103,12 @@ CCTK_REAL SpatialBetaDriverRadius "Radius at which the BetaDriver starts to be r
 } 1000000000000
 
 restricted:
+CCTK_REAL SpatialShiftGammaCoeffRadius "Radius at which the ShiftGammaCoefficient starts to be reduced"
+{
+  "*:*" :: ""
+} 1000000000000
+
+restricted:
 CCTK_INT harmonicN "d/dt alpha = - f alpha^n K  (harmonic=2, 1+log=1)"
 {
   *:* :: ""
diff --git a/ML_BSSN/src/ML_BSSN_RHS1.c b/ML_BSSN/src/ML_BSSN_RHS1.c
index 98f904a..932cdc1 100644
--- a/ML_BSSN/src/ML_BSSN_RHS1.c
+++ b/ML_BSSN/src/ML_BSSN_RHS1.c
@@ -486,17 +486,20 @@ void ML_BSSN_RHS1_Body(cGH const * restrict const cctkGH, int const dir, int con
     CCTK_REAL eta = BetaDriver*IfThen(rL > 
       SpatialBetaDriverRadius,SpatialBetaDriverRadius*INV(rL),1);
     
+    CCTK_REAL theta = IfThen(rL > SpatialShiftGammaCoeffRadius,exp(1 - 
+      rL*INV(SpatialShiftGammaCoeffRadius)),1);
+    
     CCTK_REAL beta1rhsL = (PDupwindNth1(beta1, i, j, k)*beta1L + 
       PDupwindNth2(beta1, i, j, k)*beta2L + PDupwindNth3(beta1, i, j, 
-      k)*beta3L)*ShiftAdvectionCoeff + B1L*ShiftGammaCoeff;
+      k)*beta3L)*ShiftAdvectionCoeff + B1L*ShiftGammaCoeff*theta;
     
     CCTK_REAL beta2rhsL = (PDupwindNth1(beta2, i, j, k)*beta1L + 
       PDupwindNth2(beta2, i, j, k)*beta2L + PDupwindNth3(beta2, i, j, 
-      k)*beta3L)*ShiftAdvectionCoeff + B2L*ShiftGammaCoeff;
+      k)*beta3L)*ShiftAdvectionCoeff + B2L*ShiftGammaCoeff*theta;
     
     CCTK_REAL beta3rhsL = (PDupwindNth1(beta3, i, j, k)*beta1L + 
       PDupwindNth2(beta3, i, j, k)*beta2L + PDupwindNth3(beta3, i, j, 
-      k)*beta3L)*ShiftAdvectionCoeff + B3L*ShiftGammaCoeff;
+      k)*beta3L)*ShiftAdvectionCoeff + B3L*ShiftGammaCoeff*theta;
     
     CCTK_REAL B1rhsL = -(B1L*eta) + ((PDupwindNth1(B1, i, j, k) - 
       PDupwindNth1(Xt1, i, j, k))*beta1L + (PDupwindNth2(B1, i, j, k) - 
diff --git a/ML_BSSN/src/ML_BSSN_convertFromADMBaseGamma.c b/ML_BSSN/src/ML_BSSN_convertFromADMBaseGamma.c
index 7d9e3aa..1a6159a 100644
--- a/ML_BSSN/src/ML_BSSN_convertFromADMBaseGamma.c
+++ b/ML_BSSN/src/ML_BSSN_convertFromADMBaseGamma.c
@@ -119,6 +119,7 @@ void ML_BSSN_convertFromADMBaseGamma_Body(cGH const * restrict const cctkGH, int
     CCTK_REAL  gt22L = gt22[index];
     CCTK_REAL  gt23L = gt23[index];
     CCTK_REAL  gt33L = gt33[index];
+    CCTK_REAL  rL = r[index];
     
     /* Include user supplied include files */
     
@@ -247,32 +248,44 @@ void ML_BSSN_convertFromADMBaseGamma_Body(cGH const * restrict const cctkGH, int
     CCTK_REAL AL = -(dtalpL*(-1 + 
       LapseAdvectionCoeff)*INV(harmonicF)*pow(alphaL,-harmonicN));
     
-    CCTK_REAL B1L = 6*IfThen(ShiftGammaCoeff != 
-      0,dtbetaxL*INV(ShiftGammaCoeff),0) + IfThen(ShiftGammaCoeff != 
-      0,(dtbetaxL - PDupwindNth1(betax, i, j, 
-      k)*beta1L*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff),0) + 
-      IfThen(ShiftGammaCoeff != 0,(dtbetaxL - PDupwindNth2(betax, i, j, 
-      k)*beta2L*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff),0) + 
-      IfThen(ShiftGammaCoeff != 0,(dtbetaxL - PDupwindNth3(betax, i, j, 
-      k)*beta3L*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff),0);
-    
-    CCTK_REAL B2L = 6*IfThen(ShiftGammaCoeff != 
-      0,dtbetayL*INV(ShiftGammaCoeff),0) + IfThen(ShiftGammaCoeff != 
-      0,(dtbetayL - PDupwindNth1(betay, i, j, 
-      k)*beta1L*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff),0) + 
-      IfThen(ShiftGammaCoeff != 0,(dtbetayL - PDupwindNth2(betay, i, j, 
-      k)*beta2L*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff),0) + 
-      IfThen(ShiftGammaCoeff != 0,(dtbetayL - PDupwindNth3(betay, i, j, 
-      k)*beta3L*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff),0);
-    
-    CCTK_REAL B3L = 6*IfThen(ShiftGammaCoeff != 
-      0,dtbetazL*INV(ShiftGammaCoeff),0) + IfThen(ShiftGammaCoeff != 
-      0,(dtbetazL - PDupwindNth1(betaz, i, j, 
-      k)*beta1L*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff),0) + 
-      IfThen(ShiftGammaCoeff != 0,(dtbetazL - PDupwindNth2(betaz, i, j, 
-      k)*beta2L*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff),0) + 
-      IfThen(ShiftGammaCoeff != 0,(dtbetazL - PDupwindNth3(betaz, i, j, 
-      k)*beta3L*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff),0);
+    CCTK_REAL theta = IfThen(rL > SpatialShiftGammaCoeffRadius,exp(1 - 
+      rL*INV(SpatialShiftGammaCoeffRadius)),1);
+    
+    CCTK_REAL B1L = 6*IfThen(ShiftGammaCoeff*theta != 
+      0,dtbetaxL*INV(ShiftGammaCoeff)*INV(theta),0) + 
+      IfThen(ShiftGammaCoeff*theta != 0,(dtbetaxL - PDupwindNth1(betax, i, 
+      j, 
+      k)*beta1L*ShiftAdvectionCoeff*theta)*INV(ShiftGammaCoeff)*INV(theta),0) 
+      + IfThen(ShiftGammaCoeff*theta != 0,(dtbetaxL - PDupwindNth2(betax, i, 
+      j, 
+      k)*beta2L*ShiftAdvectionCoeff*theta)*INV(ShiftGammaCoeff)*INV(theta),0) 
+      + IfThen(ShiftGammaCoeff*theta != 0,(dtbetaxL - PDupwindNth3(betax, i, 
+      j, 
+      k)*beta3L*ShiftAdvectionCoeff*theta)*INV(ShiftGammaCoeff)*INV(theta),0);
+    
+    CCTK_REAL B2L = 6*IfThen(ShiftGammaCoeff*theta != 
+      0,dtbetayL*INV(ShiftGammaCoeff)*INV(theta),0) + 
+      IfThen(ShiftGammaCoeff*theta != 0,(dtbetayL - PDupwindNth1(betay, i, 
+      j, 
+      k)*beta1L*ShiftAdvectionCoeff*theta)*INV(ShiftGammaCoeff)*INV(theta),0) 
+      + IfThen(ShiftGammaCoeff*theta != 0,(dtbetayL - PDupwindNth2(betay, i, 
+      j, 
+      k)*beta2L*ShiftAdvectionCoeff*theta)*INV(ShiftGammaCoeff)*INV(theta),0) 
+      + IfThen(ShiftGammaCoeff*theta != 0,(dtbetayL - PDupwindNth3(betay, i, 
+      j, 
+      k)*beta3L*ShiftAdvectionCoeff*theta)*INV(ShiftGammaCoeff)*INV(theta),0);
+    
+    CCTK_REAL B3L = 6*IfThen(ShiftGammaCoeff*theta != 
+      0,dtbetazL*INV(ShiftGammaCoeff)*INV(theta),0) + 
+      IfThen(ShiftGammaCoeff*theta != 0,(dtbetazL - PDupwindNth1(betaz, i, 
+      j, 
+      k)*beta1L*ShiftAdvectionCoeff*theta)*INV(ShiftGammaCoeff)*INV(theta),0) 
+      + IfThen(ShiftGammaCoeff*theta != 0,(dtbetazL - PDupwindNth2(betaz, i, 
+      j, 
+      k)*beta2L*ShiftAdvectionCoeff*theta)*INV(ShiftGammaCoeff)*INV(theta),0) 
+      + IfThen(ShiftGammaCoeff*theta != 0,(dtbetazL - PDupwindNth3(betaz, i, 
+      j, 
+      k)*beta3L*ShiftAdvectionCoeff*theta)*INV(ShiftGammaCoeff)*INV(theta),0);
     
     
     /* Copy local copies back to grid functions */
diff --git a/ML_BSSN/src/ML_BSSN_convertToADMBaseDtLapseShift.c b/ML_BSSN/src/ML_BSSN_convertToADMBaseDtLapseShift.c
index 15dbd8a..48f6554 100644
--- a/ML_BSSN/src/ML_BSSN_convertToADMBaseDtLapseShift.c
+++ b/ML_BSSN/src/ML_BSSN_convertToADMBaseDtLapseShift.c
@@ -92,6 +92,7 @@ void ML_BSSN_convertToADMBaseDtLapseShift_Body(cGH const * restrict const cctkGH
     CCTK_REAL  beta1L = beta1[index];
     CCTK_REAL  beta2L = beta2[index];
     CCTK_REAL  beta3L = beta3[index];
+    CCTK_REAL  rL = r[index];
     CCTK_REAL  trKL = trK[index];
     
     /* Include user supplied include files */
@@ -105,6 +106,9 @@ void ML_BSSN_convertToADMBaseDtLapseShift_Body(cGH const * restrict const cctkGH
     
     int dir3 = Sign(beta3L);
     
+    CCTK_REAL theta = IfThen(rL > SpatialShiftGammaCoeffRadius,exp(1 - 
+      rL*INV(SpatialShiftGammaCoeffRadius)),1);
+    
     CCTK_REAL dtalpL = (PDupwindNth1(alpha, i, j, k)*beta1L + 
       PDupwindNth2(alpha, i, j, k)*beta2L + PDupwindNth3(alpha, i, j, 
       k)*beta3L)*LapseAdvectionCoeff + harmonicF*(AL*(-1 + 
@@ -112,15 +116,15 @@ void ML_BSSN_convertToADMBaseDtLapseShift_Body(cGH const * restrict const cctkGH
     
     CCTK_REAL dtbetaxL = (PDupwindNth1(beta1, i, j, k)*beta1L + 
       PDupwindNth2(beta1, i, j, k)*beta2L + PDupwindNth3(beta1, i, j, 
-      k)*beta3L)*ShiftAdvectionCoeff + B1L*ShiftGammaCoeff;
+      k)*beta3L)*ShiftAdvectionCoeff + B1L*ShiftGammaCoeff*theta;
     
     CCTK_REAL dtbetayL = (PDupwindNth1(beta2, i, j, k)*beta1L + 
       PDupwindNth2(beta2, i, j, k)*beta2L + PDupwindNth3(beta2, i, j, 
-      k)*beta3L)*ShiftAdvectionCoeff + B2L*ShiftGammaCoeff;
+      k)*beta3L)*ShiftAdvectionCoeff + B2L*ShiftGammaCoeff*theta;
     
     CCTK_REAL dtbetazL = (PDupwindNth1(beta3, i, j, k)*beta1L + 
       PDupwindNth2(beta3, i, j, k)*beta2L + PDupwindNth3(beta3, i, j, 
-      k)*beta3L)*ShiftAdvectionCoeff + B3L*ShiftGammaCoeff;
+      k)*beta3L)*ShiftAdvectionCoeff + B3L*ShiftGammaCoeff*theta;
     
     
     /* Copy local copies back to grid functions */
diff --git a/ML_BSSN/src/ML_BSSN_convertToADMBaseDtLapseShiftBoundary.c b/ML_BSSN/src/ML_BSSN_convertToADMBaseDtLapseShiftBoundary.c
index 44874d5..40804f6 100644
--- a/ML_BSSN/src/ML_BSSN_convertToADMBaseDtLapseShiftBoundary.c
+++ b/ML_BSSN/src/ML_BSSN_convertToADMBaseDtLapseShiftBoundary.c
@@ -89,6 +89,7 @@ void ML_BSSN_convertToADMBaseDtLapseShiftBoundary_Body(cGH const * restrict cons
     CCTK_REAL  B1L = B1[index];
     CCTK_REAL  B2L = B2[index];
     CCTK_REAL  B3L = B3[index];
+    CCTK_REAL  rL = r[index];
     CCTK_REAL  trKL = trK[index];
     
     /* Include user supplied include files */
@@ -96,14 +97,17 @@ void ML_BSSN_convertToADMBaseDtLapseShiftBoundary_Body(cGH const * restrict cons
     /* Precompute derivatives */
     
     /* Calculate temporaries and grid functions */
+    CCTK_REAL theta = IfThen(rL > SpatialShiftGammaCoeffRadius,exp(1 - 
+      rL*INV(SpatialShiftGammaCoeffRadius)),1);
+    
     CCTK_REAL dtalpL = harmonicF*(AL*(-1 + LapseAdvectionCoeff) - 
       LapseAdvectionCoeff*trKL)*pow(alphaL,harmonicN);
     
-    CCTK_REAL dtbetaxL = B1L*ShiftGammaCoeff;
+    CCTK_REAL dtbetaxL = B1L*ShiftGammaCoeff*theta;
     
-    CCTK_REAL dtbetayL = B2L*ShiftGammaCoeff;
+    CCTK_REAL dtbetayL = B2L*ShiftGammaCoeff*theta;
     
-    CCTK_REAL dtbetazL = B3L*ShiftGammaCoeff;
+    CCTK_REAL dtbetazL = B3L*ShiftGammaCoeff*theta;
     
     
     /* Copy local copies back to grid functions */
diff --git a/ML_BSSN/src/ML_BSSN_convertToADMBaseFakeDtLapseShift.c b/ML_BSSN/src/ML_BSSN_convertToADMBaseFakeDtLapseShift.c
index b0f579b..828f254 100644
--- a/ML_BSSN/src/ML_BSSN_convertToADMBaseFakeDtLapseShift.c
+++ b/ML_BSSN/src/ML_BSSN_convertToADMBaseFakeDtLapseShift.c
@@ -89,6 +89,7 @@ void ML_BSSN_convertToADMBaseFakeDtLapseShift_Body(cGH const * restrict const cc
     CCTK_REAL  B1L = B1[index];
     CCTK_REAL  B2L = B2[index];
     CCTK_REAL  B3L = B3[index];
+    CCTK_REAL  rL = r[index];
     CCTK_REAL  trKL = trK[index];
     
     /* Include user supplied include files */
@@ -96,14 +97,17 @@ void ML_BSSN_convertToADMBaseFakeDtLapseShift_Body(cGH const * restrict const cc
     /* Precompute derivatives */
     
     /* Calculate temporaries and grid functions */
+    CCTK_REAL theta = IfThen(rL > SpatialShiftGammaCoeffRadius,exp(1 - 
+      rL*INV(SpatialShiftGammaCoeffRadius)),1);
+    
     CCTK_REAL dtalpL = harmonicF*(AL*(-1 + LapseAdvectionCoeff) - 
       LapseAdvectionCoeff*trKL)*pow(alphaL,harmonicN);
     
-    CCTK_REAL dtbetaxL = B1L*ShiftGammaCoeff;
+    CCTK_REAL dtbetaxL = B1L*ShiftGammaCoeff*theta;
     
-    CCTK_REAL dtbetayL = B2L*ShiftGammaCoeff;
+    CCTK_REAL dtbetayL = B2L*ShiftGammaCoeff*theta;
     
-    CCTK_REAL dtbetazL = B3L*ShiftGammaCoeff;
+    CCTK_REAL dtbetazL = B3L*ShiftGammaCoeff*theta;
     
     
     /* Copy local copies back to grid functions */