TwoPunctures: improve search for ADM masses of black holes

Incorporates the additions by Barry Wardell:

1) When the bare masses are not given as parameters (give_bare_mass=no), they
are calculated at double precision, but the par_m_plus and par_m_minus
parameters are then only set at single precision.

2) When the bare masses have been calculated, they are printed to stdout. In
some cases, the initial guess for the bare mass is printed rather than the
final calculated bare mass.

3) After looking at this section of code further, it seems that the termination
criterion for the bare mass search could be improved. Attached is an updated
patch which implements and improved check. It now continues trying to find the
bare masses until the ADM masses are within a specified tolerance of their
target. Other than the change in termination criterion, the algorithm is
unchanged, although I have simplified the calculation of a new bare mass guess
somewhat.




git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/TwoPunctures/trunk@103 b2a53a04-0f4f-0410-87ed-f9f25ced00cf

3 files changed, 74 insertions, 47 deletions

diff --git a/param.ccl b/param.ccl
index 05e5590..1978c35 100644
--- a/param.ccl
+++ b/param.ccl
@@ -39,6 +39,10 @@ BOOLEAN give_bare_mass "User provides bare masses rather than target ADM masses"
 {
 } "yes"
 
+CCTK_REAL adm_tol "Tolerance of ADM masses when give_bare_mass=no"
+{
+  (0:*) :: ""
+} 1.0e-10
 
 KEYWORD grid_setup_method "How to fill the 3D grid from the spectral grid"
 {
diff --git a/src/Newton.c b/src/Newton.c
index 54fad9e..a33b564 100644
--- a/src/Newton.c
+++ b/src/Newton.c
@@ -518,8 +518,12 @@ Newton (CCTK_POINTER_TO_CONST const cctkGH,
 #pragma omp parallel for
     for (int j = 0; j < ntotal; j++)
       dv.d0[j] = 0;
-    printf ("Newton: it=%d \t |F|=%e\n", it, (double)dmax);
-    printf ("bare mass: mp=%g \t mm=%g\n", (double) par_m_plus, (double) par_m_minus);
+
+    if(verbose==1){
+      printf ("Newton: it=%d \t |F|=%e\n", it, (double)dmax);
+      printf ("bare mass: mp=%g \t mm=%g\n", (double) par_m_plus, (double) par_m_minus);
+    }
+
     fflush(stdout);
     ii =
       bicgstab (cctkGH,
@@ -536,7 +540,10 @@ Newton (CCTK_POINTER_TO_CONST const cctkGH,
       F_of_v (cctkGH, nvar, n1, n2, n3, v, F, u);
       dmax = norm_inf (F, ntotal);
   }
-  printf ("Newton: it=%d \t |F|=%e \n", it, (double)dmax);
+
+  if(verbose==1)
+    printf ("Newton: it=%d \t |F|=%e \n", it, (double)dmax);
+
   fflush(stdout);
 
   free_dvector (F, 0, ntotal - 1);
diff --git a/src/TwoPunctures.c b/src/TwoPunctures.c
index e9284d4..e503f23 100644
--- a/src/TwoPunctures.c
+++ b/src/TwoPunctures.c
@@ -216,7 +216,7 @@ TwoPunctures (CCTK_ARGUMENTS)
 #endif
   static CCTK_REAL *F = NULL;
   static derivs u, v;
-  CCTK_REAL admMass, M_p, M_m, tmp_p, tmp_m, um, up, new_mass, old_mass;
+  CCTK_REAL admMass;
 
   if (! F) {
     /* Solve only when called for the first time */
@@ -251,59 +251,75 @@ TwoPunctures (CCTK_ARGUMENTS)
       set_initial_guess(cctkGH, v);
     }
 
+    /* If bare masses are not given, iteratively solve for them given the 
+       target ADM masses target_M_plus and target_M_minus and with initial 
+       guesses given by par_m_plus and par_m_minus. */
     if(!(give_bare_mass)) {
-      * mp  = target_M_plus;
-      * mm  = target_M_minus;
-      
-      M_p= target_M_plus;
-      M_m= target_M_minus;
-
-      new_mass = 0;
-      old_mass = 1.;
-      while (new_mass<old_mass) {
-	old_mass = (double)*mp + (double)*mm;
-	Newton (cctkGH, nvar, n1, n2, n3, v, Newton_tol, 1);
-
-	F_of_v (cctkGH, nvar, n1, n2, n3, v, F, u);
-
-	up  = PunctIntPolAtArbitPosition 
-	  (0, nvar, n1, n2, n3, v, par_b, 0., 0.);
-	um  = PunctIntPolAtArbitPosition 
-	  (0, nvar, n1, n2, n3, v, -par_b, 0., 0.);
-
-	for(int l=0; l<32;l++) {
-	  tmp_p=(4*par_b*um+ *mp +4*par_b);
-	  tmp_m=(4*par_b*up+ *mm +4*par_b);
-	  *mp = *mp - ((*mp*(up+1+M_m/tmp_p)-M_p)/
-		       (up+1+(M_m/tmp_p)-(*mp*M_m)/(pow(tmp_p,2))));
-	  *mm = *mm - ((*mm*(um+1+M_p/tmp_m)-M_m)/
-		       (um+1+(M_p/tmp_m)-(*mm*M_p)/(pow(tmp_m,2))));
-	}
-	char valbuf_p[100], valbuf_m[100];  
-	
-	sprintf (valbuf_p,"%f", (float) *mp);
-	CCTK_ParameterSet ("par_m_plus", "twopunctures", valbuf_p);
-	sprintf (valbuf_m,"%f", (float) *mm);
-	CCTK_ParameterSet ("par_m_minus", "twopunctures", valbuf_m);
-	new_mass = (double)*mp + (double)*mm;
-      }
-    } 
-    admMass = (*mp + *mm
-               - 4*par_b*PunctEvalAtArbitPosition(v.d0, 1, 0, 0, n1, n2, n3));
-    CCTK_VInfo (CCTK_THORNSTRING, "ADM mass is %g\n", (double)admMass);
+      CCTK_REAL tmp, Mp_adm, Mm_adm, Mp_adm_err, Mm_adm_err, up, um;
+      char valbuf[100];
+
+      CCTK_REAL M_p = target_M_plus;
+      CCTK_REAL M_m = target_M_minus;
+
+      CCTK_VInfo (CCTK_THORNSTRING, "Attempting to find bare masses.");
+      CCTK_VInfo (CCTK_THORNSTRING, "Target ADM masses: M_p=%g and M_m=%g",
+                  (double) M_p, (double) M_m);
+      CCTK_VInfo (CCTK_THORNSTRING, "ADM mass tolerance: %g", (double) adm_tol);
+
+      /* Loop until both ADM masses are within adm_tol of their target */
+      do {
+        CCTK_VInfo (CCTK_THORNSTRING, "Bare masses: mp=%g, mm=%g",
+                    (double)*mp, (double)*mm);
+        Newton (cctkGH, nvar, n1, n2, n3, v, Newton_tol, 1);
+
+        F_of_v (cctkGH, nvar, n1, n2, n3, v, F, u);
+
+        up = PunctIntPolAtArbitPosition(0, nvar, n1, n2, n3, v, par_b, 0., 0.);
+        um = PunctIntPolAtArbitPosition(0, nvar, n1, n2, n3, v,-par_b, 0., 0.);
+
+        /* Calculate the ADM masses from the current bare mass guess */
+        Mp_adm = (1 + up) * *mp + *mp * *mm / (4. * par_b);
+        Mm_adm = (1 + um) * *mm + *mp * *mm / (4. * par_b);
+
+        /* Check how far the current ADM masses are from the target */
+        Mp_adm_err = fabs(M_p-Mp_adm);
+        Mm_adm_err = fabs(M_m-Mm_adm);
+        CCTK_VInfo (CCTK_THORNSTRING, "ADM mass error: M_p_err=%.4g, M_m_err=%.4g",
+                    (double)Mp_adm_err, (double)Mm_adm_err);
+        
+        /* Invert the ADM mass equation and update the bare mass guess so that
+           it gives the correct target ADM masses */
+        tmp = -4*par_b*( 1 + um + up + um*up ) + 
+                sqrt(16*par_b*M_m*(1 + um)*(1 + up) + 
+                  pow(-M_m + M_p + 4*par_b*(1 + um)*(1 + up),2));
+        *mp = (tmp + M_p - M_m)/(2.*(1 + up));
+        *mm = (tmp - M_p + M_m)/(2.*(1 + um));
+        
+        /* Set the par_m_plus and par_m_minus parameters */
+        sprintf (valbuf,"%.17g", (double) *mp);
+        CCTK_ParameterSet ("par_m_plus", "twopunctures", valbuf);
+        
+        sprintf (valbuf,"%.17g", (double) *mm);
+        CCTK_ParameterSet ("par_m_minus", "twopunctures", valbuf);
+        
+      } while ( (Mp_adm_err > adm_tol) ||
+                (Mm_adm_err > adm_tol) );
+                
+      CCTK_VInfo (CCTK_THORNSTRING, "Found bare masses.");
+    }
 
     Newton (cctkGH, nvar, n1, n2, n3, v, Newton_tol, Newton_maxit);
  
     F_of_v (cctkGH, nvar, n1, n2, n3, v, F, u);
 
     CCTK_VInfo (CCTK_THORNSTRING,
-		  "The two puncture masses are %g and %g",
-                (double) par_m_minus, (double) par_m_plus);
+		  "The two puncture masses are mp=%.17g and mm=%.17g",
+                (double) *mp, (double) *mm);
 
     /* print out ADM mass, eq.: \Delta M_ADM=2*r*u=4*b*V for A=1,B=0,phi=0 */
     admMass = (*mp + *mm
-               - 4*par_b*PunctEvalAtArbitPosition(v.d0, 1, 0, 0, n1, n2, n3));
-    CCTK_VInfo (CCTK_THORNSTRING, "ADM mass is %g", (double)admMass);
+               - 4*par_b*PunctEvalAtArbitPosition(v.d0, 1, 0, 0, n1, n2, n3));;
+    CCTK_VInfo (CCTK_THORNSTRING, "The total ADM mass is %g", (double) admMass);
   }
 
   if (CCTK_EQUALS(grid_setup_method, "Taylor expansion"))