Use a fudge factor to remove singularities when grid points fall onto

an axis.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/IDAnalyticBH/trunk@152 6a3ddf76-46e1-4315-99d9-bc56cac1ef84

2 files changed, 19 insertions, 8 deletions

diff --git a/param.ccl b/param.ccl
index 933c261..303531e 100644
--- a/param.ccl
+++ b/param.ccl
@@ -17,7 +17,6 @@ REAL a_Kerr "Angular momentum parameter of black hole"
 -1:1 :: "Between +1 and -1"
 } 0.1
 
-
 # Multiple Misner Parameters
 # --------------------------
 REAL mu "Misner mu value"
@@ -106,6 +105,13 @@ REAL bl_M_4 "Mass of 4th BL hole"
 : :: "Anything"
 } 1.0
 
+# Common parameters
+# -----------------
+REAL epsilon "Fudge factor"
+{
+  0.0:* :: ""
+} 1.e-16
+
 
 shares: ADMBase
 
diff --git a/src/Schwarzschild.c b/src/Schwarzschild.c
index 850bd9a..331b10a 100644
--- a/src/Schwarzschild.c
+++ b/src/Schwarzschild.c
@@ -9,6 +9,7 @@
    @version   $Id$
  @@*/
 
+#include <math.h>
 #include <string.h>
 
 #include "cctk.h"
@@ -68,14 +69,16 @@ void Schwarzschild(CCTK_ARGUMENTS)
     for (i = 0; i < npoints; i++)
     {
       /*        Compute conformal factor */
-      psi[i] = ( one + mass/(two*r[i]+1e-20));
+      CCTK_REAL const rr = pow(pow(r[i], 4) + pow(epsilon, 4), 0.25);
+
+      psi[i] = (one + mass / (two*rr));
 
       if(make_conformal_derivs)
       {
         /*        derivatives of psi / psi */
-        r_squared = r[i]*r[i];
-        r_cubed   = r[i]*r_squared;
-        tmp = mass/(two*r_cubed*psi[i]+1e-20);
+        r_squared = rr*rr;
+        r_cubed   = rr*r_squared;
+        tmp = mass / (two*r_cubed*psi[i]);
 
 
         psix[i] = -x[i]*tmp;
@@ -84,7 +87,7 @@ void Schwarzschild(CCTK_ARGUMENTS)
 
         if(*conformal_state > 2)
         {
-          tmp = mass/(two*r_squared*r_cubed*psi[i]+1e-20);
+          tmp = mass/(two*r_squared*r_cubed*psi[i]);
           psixy[i] = three*x[i]*y[i]*tmp;
           psixz[i] = three*x[i]*z[i]*tmp;
           psiyz[i] = three*y[i]*z[i]*tmp;
@@ -106,7 +109,9 @@ void Schwarzschild(CCTK_ARGUMENTS)
   {
     for (i = 0; i < npoints; i++)
     {
-      tmp = one + mass/(two*r[i]+1e-20);
+      CCTK_REAL const rr = pow(pow(r[i], 4) + pow(epsilon, 4), 0.25);
+
+      tmp = one + mass / (two*rr);
       gxx[i] = tmp*tmp*tmp*tmp;
       gyy[i] = gxx[i];
       gzz[i] = gxx[i];
@@ -124,7 +129,7 @@ void Schwarzschild(CCTK_ARGUMENTS)
 
     for (i = 0; i < npoints; i++)
     {
-      alp[i] = (2.*r[i] - mass)/(2.*r[i]+mass);
+      alp[i] = (two*r[i] - mass) / (two*r[i] + mass);
     }
   }