Avoid division-by-zero problems when ShiftGammaCoeff=0

7 files changed, 93 insertions, 57 deletions

diff --git a/ML_BSSN/src/ML_BSSN_convertFromADMBaseGamma.c b/ML_BSSN/src/ML_BSSN_convertFromADMBaseGamma.c
index ef2a0a3..8de9e36 100644
--- a/ML_BSSN/src/ML_BSSN_convertFromADMBaseGamma.c
+++ b/ML_BSSN/src/ML_BSSN_convertFromADMBaseGamma.c
@@ -266,14 +266,23 @@ void ML_BSSN_convertFromADMBaseGamma_Body(cGH const * const cctkGH, CCTK_INT con
     
     AL  =  -(dtalpL*(-1 + LapseAdvectionCoeff)*INV(harmonicF)*pow(alphaL,-harmonicN));
     
-    B1L  =  (dtbetaxL - (PDupwindNth1(beta1, i, j, k)*beta1L + PDupwindNth2(beta1, i, j, k)*beta2L + 
-             PDupwindNth3(beta1, i, j, k)*beta3L)*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff);
-    
-    B2L  =  (dtbetayL - (PDupwindNth1(beta2, i, j, k)*beta1L + PDupwindNth2(beta2, i, j, k)*beta2L + 
-             PDupwindNth3(beta2, i, j, k)*beta3L)*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff);
-    
-    B3L  =  (dtbetazL - (PDupwindNth1(beta3, i, j, k)*beta1L + PDupwindNth2(beta3, i, j, k)*beta2L + 
-             PDupwindNth3(beta3, i, j, k)*beta3L)*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff);
+    B1L  =  6*IfThen(ShiftGammaCoeff,dtbetaxL*INV(ShiftGammaCoeff),0) + 
+        IfThen(ShiftGammaCoeff,(dtbetaxL - PDupwindNth1(beta1, i, j, k)*beta1L*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff),
+         0) + IfThen(ShiftGammaCoeff,(dtbetaxL - PDupwindNth2(beta1, i, j, k)*beta2L*ShiftAdvectionCoeff)*
+          INV(ShiftGammaCoeff),0) + IfThen(ShiftGammaCoeff,
+         (dtbetaxL - PDupwindNth3(beta1, i, j, k)*beta3L*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff),0);
+    
+    B2L  =  6*IfThen(ShiftGammaCoeff,dtbetayL*INV(ShiftGammaCoeff),0) + 
+        IfThen(ShiftGammaCoeff,(dtbetayL - PDupwindNth1(beta2, i, j, k)*beta1L*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff),
+         0) + IfThen(ShiftGammaCoeff,(dtbetayL - PDupwindNth2(beta2, i, j, k)*beta2L*ShiftAdvectionCoeff)*
+          INV(ShiftGammaCoeff),0) + IfThen(ShiftGammaCoeff,
+         (dtbetayL - PDupwindNth3(beta2, i, j, k)*beta3L*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff),0);
+    
+    B3L  =  6*IfThen(ShiftGammaCoeff,dtbetazL*INV(ShiftGammaCoeff),0) + 
+        IfThen(ShiftGammaCoeff,(dtbetazL - PDupwindNth1(beta3, i, j, k)*beta1L*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff),
+         0) + IfThen(ShiftGammaCoeff,(dtbetazL - PDupwindNth2(beta3, i, j, k)*beta2L*ShiftAdvectionCoeff)*
+          INV(ShiftGammaCoeff),0) + IfThen(ShiftGammaCoeff,
+         (dtbetazL - PDupwindNth3(beta3, i, j, k)*beta3L*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff),0);
     
     
     /* Copy local copies back to grid functions */
diff --git a/ML_BSSN_Helper/schedule.ccl b/ML_BSSN_Helper/schedule.ccl
index 447b2cd..97c5c8c 100644
--- a/ML_BSSN_Helper/schedule.ccl
+++ b/ML_BSSN_Helper/schedule.ccl
@@ -118,13 +118,6 @@ if (CCTK_EQUALS (evolution_method, "ML_BSSN")) {
   
   if (CCTK_EQUALS (calculate_ADMBase_variables_at, "MoL_PostStep"))
   {
-    # if (timelevels > 1)
-    # {
-    #   SCHEDULE ML_BSSN_CopyADMBase AT prestep
-    #   {
-    #     LANG: C
-    #   } "Copy ADMBase variables to current time level"
-    # }
     SCHEDULE GROUP ML_BSSN_convertToADMBaseGroupWrapper IN MoL_PostStep AFTER (ML_BSSN_ApplyBCs ML_BSSN_enforce) BEFORE (ADMBase_SetADMVars Whisky_PostStep)
     {
     } "Calculate ADM variables"
diff --git a/ML_BSSN_M/src/ML_BSSN_M_convertFromADMBaseGamma.c b/ML_BSSN_M/src/ML_BSSN_M_convertFromADMBaseGamma.c
index 3e2f266..1b93415 100644
--- a/ML_BSSN_M/src/ML_BSSN_M_convertFromADMBaseGamma.c
+++ b/ML_BSSN_M/src/ML_BSSN_M_convertFromADMBaseGamma.c
@@ -266,14 +266,23 @@ void ML_BSSN_M_convertFromADMBaseGamma_Body(cGH const * const cctkGH, CCTK_INT c
     
     AL  =  -(dtalpL*(-1 + LapseAdvectionCoeff)*INV(harmonicF)*pow(alphaL,-harmonicN));
     
-    B1L  =  (dtbetaxL - (PDupwindNth1(beta1, i, j, k)*beta1L + PDupwindNth2(beta1, i, j, k)*beta2L + 
-             PDupwindNth3(beta1, i, j, k)*beta3L)*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff);
-    
-    B2L  =  (dtbetayL - (PDupwindNth1(beta2, i, j, k)*beta1L + PDupwindNth2(beta2, i, j, k)*beta2L + 
-             PDupwindNth3(beta2, i, j, k)*beta3L)*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff);
-    
-    B3L  =  (dtbetazL - (PDupwindNth1(beta3, i, j, k)*beta1L + PDupwindNth2(beta3, i, j, k)*beta2L + 
-             PDupwindNth3(beta3, i, j, k)*beta3L)*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff);
+    B1L  =  6*IfThen(ShiftGammaCoeff,dtbetaxL*INV(ShiftGammaCoeff),0) + 
+        IfThen(ShiftGammaCoeff,(dtbetaxL - PDupwindNth1(beta1, i, j, k)*beta1L*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff),
+         0) + IfThen(ShiftGammaCoeff,(dtbetaxL - PDupwindNth2(beta1, i, j, k)*beta2L*ShiftAdvectionCoeff)*
+          INV(ShiftGammaCoeff),0) + IfThen(ShiftGammaCoeff,
+         (dtbetaxL - PDupwindNth3(beta1, i, j, k)*beta3L*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff),0);
+    
+    B2L  =  6*IfThen(ShiftGammaCoeff,dtbetayL*INV(ShiftGammaCoeff),0) + 
+        IfThen(ShiftGammaCoeff,(dtbetayL - PDupwindNth1(beta2, i, j, k)*beta1L*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff),
+         0) + IfThen(ShiftGammaCoeff,(dtbetayL - PDupwindNth2(beta2, i, j, k)*beta2L*ShiftAdvectionCoeff)*
+          INV(ShiftGammaCoeff),0) + IfThen(ShiftGammaCoeff,
+         (dtbetayL - PDupwindNth3(beta2, i, j, k)*beta3L*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff),0);
+    
+    B3L  =  6*IfThen(ShiftGammaCoeff,dtbetazL*INV(ShiftGammaCoeff),0) + 
+        IfThen(ShiftGammaCoeff,(dtbetazL - PDupwindNth1(beta3, i, j, k)*beta1L*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff),
+         0) + IfThen(ShiftGammaCoeff,(dtbetazL - PDupwindNth2(beta3, i, j, k)*beta2L*ShiftAdvectionCoeff)*
+          INV(ShiftGammaCoeff),0) + IfThen(ShiftGammaCoeff,
+         (dtbetazL - PDupwindNth3(beta3, i, j, k)*beta3L*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff),0);
     
     
     /* Copy local copies back to grid functions */
diff --git a/ML_BSSN_MP/src/ML_BSSN_MP_convertFromADMBaseGamma.c b/ML_BSSN_MP/src/ML_BSSN_MP_convertFromADMBaseGamma.c
index a7c4713..b89269c 100644
--- a/ML_BSSN_MP/src/ML_BSSN_MP_convertFromADMBaseGamma.c
+++ b/ML_BSSN_MP/src/ML_BSSN_MP_convertFromADMBaseGamma.c
@@ -328,20 +328,50 @@ void ML_BSSN_MP_convertFromADMBaseGamma_Body(cGH const * const cctkGH, CCTK_INT
     
     AL  =  -(dtalpL*(-1 + LapseAdvectionCoeff)*INV(harmonicF)*pow(alphaL,-harmonicN));
     
-    B1L  =  (dtbetaxL - (PDupwindNth1(beta1, i, j, k)*(beta1L*J11L + beta2L*J12L + beta3L*J13L) + 
-             PDupwindNth2(beta1, i, j, k)*(beta1L*J21L + beta2L*J22L + beta3L*J23L) + 
-             PDupwindNth3(beta1, i, j, k)*(beta1L*J31L + beta2L*J32L + beta3L*J33L))*ShiftAdvectionCoeff)*
-        INV(ShiftGammaCoeff);
-    
-    B2L  =  (dtbetayL - (PDupwindNth1(beta2, i, j, k)*(beta1L*J11L + beta2L*J12L + beta3L*J13L) + 
-             PDupwindNth2(beta2, i, j, k)*(beta1L*J21L + beta2L*J22L + beta3L*J23L) + 
-             PDupwindNth3(beta2, i, j, k)*(beta1L*J31L + beta2L*J32L + beta3L*J33L))*ShiftAdvectionCoeff)*
-        INV(ShiftGammaCoeff);
-    
-    B3L  =  (dtbetazL - (PDupwindNth1(beta3, i, j, k)*(beta1L*J11L + beta2L*J12L + beta3L*J13L) + 
-             PDupwindNth2(beta3, i, j, k)*(beta1L*J21L + beta2L*J22L + beta3L*J23L) + 
-             PDupwindNth3(beta3, i, j, k)*(beta1L*J31L + beta2L*J32L + beta3L*J33L))*ShiftAdvectionCoeff)*
-        INV(ShiftGammaCoeff);
+    B1L  =  IfThen(ShiftGammaCoeff,(dtbetaxL - PDupwindNth1(beta1, i, j, k)*beta1L*J11L*ShiftAdvectionCoeff)*
+          INV(ShiftGammaCoeff),0) + IfThen(ShiftGammaCoeff,
+         (dtbetaxL - PDupwindNth1(beta1, i, j, k)*beta2L*J12L*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff),0) + 
+        IfThen(ShiftGammaCoeff,(dtbetaxL - PDupwindNth1(beta1, i, j, k)*beta3L*J13L*ShiftAdvectionCoeff)*
+          INV(ShiftGammaCoeff),0) + IfThen(ShiftGammaCoeff,
+         (dtbetaxL - PDupwindNth2(beta1, i, j, k)*beta1L*J21L*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff),0) + 
+        IfThen(ShiftGammaCoeff,(dtbetaxL - PDupwindNth2(beta1, i, j, k)*beta2L*J22L*ShiftAdvectionCoeff)*
+          INV(ShiftGammaCoeff),0) + IfThen(ShiftGammaCoeff,
+         (dtbetaxL - PDupwindNth2(beta1, i, j, k)*beta3L*J23L*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff),0) + 
+        IfThen(ShiftGammaCoeff,(dtbetaxL - PDupwindNth3(beta1, i, j, k)*beta1L*J31L*ShiftAdvectionCoeff)*
+          INV(ShiftGammaCoeff),0) + IfThen(ShiftGammaCoeff,
+         (dtbetaxL - PDupwindNth3(beta1, i, j, k)*beta2L*J32L*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff),0) + 
+        IfThen(ShiftGammaCoeff,(dtbetaxL - PDupwindNth3(beta1, i, j, k)*beta3L*J33L*ShiftAdvectionCoeff)*
+          INV(ShiftGammaCoeff),0);
+    
+    B2L  =  IfThen(ShiftGammaCoeff,(dtbetayL - PDupwindNth1(beta2, i, j, k)*beta1L*J11L*ShiftAdvectionCoeff)*
+          INV(ShiftGammaCoeff),0) + IfThen(ShiftGammaCoeff,
+         (dtbetayL - PDupwindNth1(beta2, i, j, k)*beta2L*J12L*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff),0) + 
+        IfThen(ShiftGammaCoeff,(dtbetayL - PDupwindNth1(beta2, i, j, k)*beta3L*J13L*ShiftAdvectionCoeff)*
+          INV(ShiftGammaCoeff),0) + IfThen(ShiftGammaCoeff,
+         (dtbetayL - PDupwindNth2(beta2, i, j, k)*beta1L*J21L*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff),0) + 
+        IfThen(ShiftGammaCoeff,(dtbetayL - PDupwindNth2(beta2, i, j, k)*beta2L*J22L*ShiftAdvectionCoeff)*
+          INV(ShiftGammaCoeff),0) + IfThen(ShiftGammaCoeff,
+         (dtbetayL - PDupwindNth2(beta2, i, j, k)*beta3L*J23L*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff),0) + 
+        IfThen(ShiftGammaCoeff,(dtbetayL - PDupwindNth3(beta2, i, j, k)*beta1L*J31L*ShiftAdvectionCoeff)*
+          INV(ShiftGammaCoeff),0) + IfThen(ShiftGammaCoeff,
+         (dtbetayL - PDupwindNth3(beta2, i, j, k)*beta2L*J32L*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff),0) + 
+        IfThen(ShiftGammaCoeff,(dtbetayL - PDupwindNth3(beta2, i, j, k)*beta3L*J33L*ShiftAdvectionCoeff)*
+          INV(ShiftGammaCoeff),0);
+    
+    B3L  =  IfThen(ShiftGammaCoeff,(dtbetazL - PDupwindNth1(beta3, i, j, k)*beta1L*J11L*ShiftAdvectionCoeff)*
+          INV(ShiftGammaCoeff),0) + IfThen(ShiftGammaCoeff,
+         (dtbetazL - PDupwindNth1(beta3, i, j, k)*beta2L*J12L*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff),0) + 
+        IfThen(ShiftGammaCoeff,(dtbetazL - PDupwindNth1(beta3, i, j, k)*beta3L*J13L*ShiftAdvectionCoeff)*
+          INV(ShiftGammaCoeff),0) + IfThen(ShiftGammaCoeff,
+         (dtbetazL - PDupwindNth2(beta3, i, j, k)*beta1L*J21L*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff),0) + 
+        IfThen(ShiftGammaCoeff,(dtbetazL - PDupwindNth2(beta3, i, j, k)*beta2L*J22L*ShiftAdvectionCoeff)*
+          INV(ShiftGammaCoeff),0) + IfThen(ShiftGammaCoeff,
+         (dtbetazL - PDupwindNth2(beta3, i, j, k)*beta3L*J23L*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff),0) + 
+        IfThen(ShiftGammaCoeff,(dtbetazL - PDupwindNth3(beta3, i, j, k)*beta1L*J31L*ShiftAdvectionCoeff)*
+          INV(ShiftGammaCoeff),0) + IfThen(ShiftGammaCoeff,
+         (dtbetazL - PDupwindNth3(beta3, i, j, k)*beta2L*J32L*ShiftAdvectionCoeff)*INV(ShiftGammaCoeff),0) + 
+        IfThen(ShiftGammaCoeff,(dtbetazL - PDupwindNth3(beta3, i, j, k)*beta3L*J33L*ShiftAdvectionCoeff)*
+          INV(ShiftGammaCoeff),0);
     
     
     /* Copy local copies back to grid functions */
diff --git a/ML_BSSN_MP_Helper/schedule.ccl b/ML_BSSN_MP_Helper/schedule.ccl
index 0d44d8d..fde3844 100644
--- a/ML_BSSN_MP_Helper/schedule.ccl
+++ b/ML_BSSN_MP_Helper/schedule.ccl
@@ -118,13 +118,6 @@ if (CCTK_EQUALS (evolution_method, "ML_BSSN_MP")) {
   
   if (CCTK_EQUALS (calculate_ADMBase_variables_at, "MoL_PostStep"))
   {
-    # if (timelevels > 1)
-    # {
-    #   SCHEDULE ML_BSSN_MP_CopyADMBase AT prestep
-    #   {
-    #     LANG: C
-    #   } "Copy ADMBase variables to current time level"
-    # }
     SCHEDULE GROUP ML_BSSN_MP_convertToADMBaseGroupWrapper IN MoL_PostStep AFTER (ML_BSSN_MP_ApplyBCs ML_BSSN_MP_enforce) BEFORE (ADMBase_SetADMVars Whisky_PostStep)
     {
     } "Calculate ADM variables"
diff --git a/ML_BSSN_M_Helper/schedule.ccl b/ML_BSSN_M_Helper/schedule.ccl
index ee7adac..ba739c5 100644
--- a/ML_BSSN_M_Helper/schedule.ccl
+++ b/ML_BSSN_M_Helper/schedule.ccl
@@ -118,13 +118,6 @@ if (CCTK_EQUALS (evolution_method, "ML_BSSN_M")) {
   
   if (CCTK_EQUALS (calculate_ADMBase_variables_at, "MoL_PostStep"))
   {
-    # if (timelevels > 1)
-    # {
-    #   SCHEDULE ML_BSSN_M_CopyADMBase AT prestep
-    #   {
-    #     LANG: C
-    #   } "Copy ADMBase variables to current time level"
-    # }
     SCHEDULE GROUP ML_BSSN_M_convertToADMBaseGroupWrapper IN MoL_PostStep AFTER (ML_BSSN_M_ApplyBCs ML_BSSN_M_enforce) BEFORE (ADMBase_SetADMVars Whisky_PostStep)
     {
     } "Calculate ADM variables"
diff --git a/m/McLachlan_BSSN.m b/m/McLachlan_BSSN.m
index da0bd80..c6a4dae 100644
--- a/m/McLachlan_BSSN.m
+++ b/m/McLachlan_BSSN.m
@@ -345,12 +345,21 @@ convertFromADMBaseGammaCalc =
     
     A -> - dtalp / (harmonicF alpha^harmonicN) (LapseAdvectionCoeff - 1),
     
-    B1 -> 1/ShiftGammaCoeff
-          (dtbetax - ShiftAdvectionCoeff beta[ua] PDu[beta1,la]),
-    B2 -> 1/ShiftGammaCoeff
-          (dtbetay - ShiftAdvectionCoeff beta[ua] PDu[beta2,la]),
-    B3 -> 1/ShiftGammaCoeff
-          (dtbetaz - ShiftAdvectionCoeff beta[ua] PDu[beta3,la])
+    (* If ShiftGammaCoeff=0, then B^i is not evolved, in the sense
+       that it does not influence the time evolution of other
+       variables.  *)
+    B1 -> IfThen[ShiftGammaCoeff,
+                 1/ShiftGammaCoeff
+                 (dtbetax - ShiftAdvectionCoeff beta[ua] PDu[beta1,la]),
+                 0],
+    B2 -> IfThen[ShiftGammaCoeff,
+                 1/ShiftGammaCoeff
+                 (dtbetay - ShiftAdvectionCoeff beta[ua] PDu[beta2,la]),
+                 0],
+    B3 -> IfThen[ShiftGammaCoeff,
+                 1/ShiftGammaCoeff
+                 (dtbetaz - ShiftAdvectionCoeff beta[ua] PDu[beta3,la]),
+                 0]
   }
 };