From 087a8ee26d8abfccaf428cb43db81a64e00efe8b Mon Sep 17 00:00:00 2001
From: knarf <knarf@b2a53a04-0f4f-0410-87ed-f9f25ced00cf>
Date: Fri, 2 Dec 2005 10:26:08 +0000
Subject: a few smaller fixes allow an extended region blended instead of the
 puncture % done - info (locally, but who cares ...)

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/TwoPunctures/trunk@48 b2a53a04-0f4f-0410-87ed-f9f25ced00cf
---
 param.ccl          |  7 +++++--
 src/TwoPunctures.c | 53 ++++++++++++++++++++++++++++++++++++++++++++++++-----
 2 files changed, 53 insertions(+), 7 deletions(-)

diff --git a/param.ccl b/param.ccl
index 22a73d9..8764a16 100644
--- a/param.ccl
+++ b/param.ccl
@@ -70,12 +70,15 @@ INT Newton_maxit "Maximum number of Newton iterations"
 {
   0:* :: ""
 } 5
+
 REAL TP_Tiny "Tiny number to avoid nans near or at the pucture locations"
 {
   0:* :: "anything positive, usually very small"
 } 0.0
-
-
+REAL TP_Extend_Radius "Radius of an extended spacetime instead of the puncture"
+{
+  0:* :: "anything positive, should be smaller than the horizon"
+} 0.0
 
 REAL par_b "x coordinate of the m+ puncture"
 {
diff --git a/src/TwoPunctures.c b/src/TwoPunctures.c
index 4ae122f..b41b135 100644
--- a/src/TwoPunctures.c
+++ b/src/TwoPunctures.c
@@ -196,6 +196,7 @@ TwoPunctures (CCTK_ARGUMENTS)
   int nvar = 1, n1 = npoints_A, n2 = npoints_B, n3 = npoints_phi;
 
   int i, j, k, ntotal = n1 * n2 * n3 * nvar;
+  int percent;
   static double *F = NULL;
   static derivs u, v;
   double admMass;
@@ -272,6 +273,13 @@ TwoPunctures (CCTK_ARGUMENTS)
     {
       for (i = 0; i < cctk_lsh[0]; ++i)
       {
+        if (percent != 10*(i+j*cctk_lsh[0]+k*cctk_lsh[0]*cctk_lsh[1]) /
+                          (cctk_lsh[0] * cctk_lsh[1] * cctk_lsh[2]))
+        {
+            percent = 10*(i+j*cctk_lsh[0]+k*cctk_lsh[0]*cctk_lsh[1]) /
+                         (cctk_lsh[0] * cctk_lsh[1] * cctk_lsh[2]);
+            CCTK_VInfo(CCTK_THORNSTRING, "%3d%% done", percent*10);
+        }
 
 	const int ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
         
@@ -298,9 +306,23 @@ TwoPunctures (CCTK_ARGUMENTS)
             r_plus = TP_Tiny;
         if (r_minus < TP_Tiny)
             r_minus = TP_Tiny;
-        const double psi1 = 1
+        double psi1 = 1
           + 0.5 * par_m_plus / r_plus
           + 0.5 * par_m_minus / r_minus + U;
+#define EXTEND(M,r) \
+          ( M * (3./8 * pow(r, 4.0) / pow(TP_Extend_Radius, 5.0) - \
+                 5./4 * pow(r, 2.0) / pow(TP_Extend_Radius, 3.0) + \
+                 15./8 / TP_Extend_Radius))
+        if (r_plus < TP_Extend_Radius) {
+          psi1 = 1
+             + 0.5 * EXTEND(par_m_plus,r_plus)
+             + 0.5 * par_m_minus / r_minus + U;
+        }
+        if (r_minus < TP_Extend_Radius) {
+          psi1 = 1
+             + 0.5 * EXTEND(par_m_minus,r_minus)
+             + 0.5 * par_m_plus / r_plus + U;
+        }
         double static_psi = 1;
         
         double Aij[3][3];
@@ -311,7 +333,6 @@ TwoPunctures (CCTK_ARGUMENTS)
           double xp, yp, zp, rp, ir;
           double s1, s3, s5;
           double p, px, py, pz, pxx, pxy, pxz, pyy, pyz, pzz;
-
           p = 1.0;
           px = py = pz = 0.0;
           pxx = pxy = pxz = 0.0;
@@ -322,8 +343,14 @@ TwoPunctures (CCTK_ARGUMENTS)
           yp = y[ind];
           zp = z[ind];
           rp = sqrt (xp*xp + yp*yp + zp*zp);
+          if (rp < TP_Tiny)
+              rp = TP_Tiny;
           ir = 1.0/rp;
 
+          if (rp < TP_Extend_Radius) {
+            ir = EXTEND(1., rp);
+          }
+
           s1 = 0.5*par_m_plus*ir;
           s3 = -s1*ir*ir;
           s5 = -3.0*s3*ir*ir;
@@ -346,8 +373,14 @@ TwoPunctures (CCTK_ARGUMENTS)
           yp = y[ind];
           zp = z[ind];
           rp = sqrt (xp*xp + yp*yp + zp*zp);
+          if (rp < TP_Tiny)
+              rp = TP_Tiny;
           ir = 1.0/rp;
 
+          if (rp < TP_Extend_Radius) {
+            ir = EXTEND(1., rp);
+          }
+
           s1 = 0.5*par_m_minus*ir;
           s3 = -s1*ir*ir;
           s5 = -3.0*s3*ir*ir;
@@ -408,14 +441,24 @@ TwoPunctures (CCTK_ARGUMENTS)
 /*                                / (1.0 + 0.5 * par_m_minus / r_minus)); */
 /*           alp[ind] = alp1 * alp2; */
           alp[ind] =
-            ((1.0 - 0.5 * par_m_plus / r_plus - 0.5 * par_m_minus / r_minus)
-             / (1.0 + 0.5 * par_m_plus / r_plus + 0.5 * par_m_minus / r_minus));
+            ((1.0 -0.5*par_m_plus/r_plus -0.5*par_m_minus/r_minus)
+            /(1.0 +0.5*par_m_plus/r_plus +0.5*par_m_minus/r_minus));
+
+          if (r_plus < TP_Extend_Radius) {
+            alp[ind] =
+              ((1.0 -0.5*EXTEND(par_m_plus, r_plus) -0.5*par_m_minus/r_minus)
+              /(1.0 +0.5*EXTEND(par_m_plus, r_plus) +0.5*par_m_minus/r_minus));
+          }
+          if (r_minus < TP_Extend_Radius) {
+            alp[ind] =
+              ((1.0 -0.5*EXTEND(par_m_minus, r_minus) -0.5*par_m_plus/r_plus)
+              /(1.0 +0.5*EXTEND(par_m_minus, r_minus) +0.5*par_m_plus/r_plus));
+          }
           
           if (averaged_lapse) {
             alp[ind] = 0.5 * (1.0 + alp[ind]);
           }
         }
-
       }
     }
   }
-- 
cgit v1.2.3