Handle the case correctly where ShiftGammaCoeff=0. Don't dividy by ShiftGammaCoeff in this case; instead, initialise the (unused) B^i to zero.

4 files changed, 95 insertions, 38 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_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/m/McLachlan_BSSN.m b/m/McLachlan_BSSN.m
index c18cd88..ba7c541 100644
--- a/m/McLachlan_BSSN.m
+++ b/m/McLachlan_BSSN.m
@@ -1,5 +1,5 @@
-$Path = Join[$Path, {"~/Calpha/kranc/Tools/CodeGen",
-                     "~/Calpha/kranc/Tools/MathematicaMisc"}];
+$Path = Join[$Path, {"../../../kranc/Tools/CodeGen",
+                     "../../../kranc/Tools/MathematicaMisc"}];
 
 Get["KrancThorn`"];
 
@@ -351,12 +351,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]
   }
 };