implement UIUC's speed-up in "evaluation" of spectral solution

The accompanying svn patch applies the recently released refinement of
TwoPunctures by Vasileios Paschalidis and Zach Etienne at UIUC, greatly
reducing the time taken to properly apply the solution of the spectral solve to
all grid points.

From the abstract of the notes in arXiv:1304.0457:

TwoPunctures is perhaps the most widely-adopted code for generating binary
black hole "puncture" initial data and interpolating these (spectral) data onto
evolution grids. In typical usage, the bulk of this code's run time is spent in
its spectral interpolation routine. We announce a new publicly-available
spectral interpolation routine that improves the performance of the original
interpolation routine by a factor of ~100, yielding results consistent with the
original spectral interpolation routine to roundoff precision. This note serves
as a guide for installing this routine both in the original standalone
TwoPunctures code and the Einstein Toolkit supported version of this code. 


Patch kindly provided by Bernard Kelly
Original code by Vasileios Paschalidis and Zach Etienne


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

1 files changed, 21 insertions, 8 deletions

diff --git a/src/TwoPunctures.c b/src/TwoPunctures.c
index 8835c74..bad5465 100644
--- a/src/TwoPunctures.c
+++ b/src/TwoPunctures.c
@@ -216,7 +216,7 @@ TwoPunctures (CCTK_ARGUMENTS)
   int percent10 = 0;
 #endif
   static CCTK_REAL *F = NULL;
-  static derivs u, v;
+  static derivs u, v, cf_v;
   CCTK_REAL admMass;
 
   if (! F) {
@@ -225,6 +225,7 @@ TwoPunctures (CCTK_ARGUMENTS)
     F = dvector (0, ntotal - 1);
     allocate_derivs (&u, ntotal);
     allocate_derivs (&v, ntotal);
+    allocate_derivs (&cf_v, ntotal);
 
     if (use_sources) {
       CCTK_INFO ("Solving puncture equation for BH-NS/NS-NS system");
@@ -236,6 +237,16 @@ TwoPunctures (CCTK_ARGUMENTS)
     /* initialise to 0 */
     for (int j = 0; j < ntotal; j++)
     {
+      cf_v.d0[j] = 0.0;
+      cf_v.d1[j] = 0.0;
+      cf_v.d2[j] = 0.0;
+      cf_v.d3[j] = 0.0;
+      cf_v.d11[j] = 0.0;
+      cf_v.d12[j] = 0.0;
+      cf_v.d13[j] = 0.0;
+      cf_v.d22[j] = 0.0;
+      cf_v.d23[j] = 0.0;
+      cf_v.d33[j] = 0.0;
       v.d0[j] = 0.0;
       v.d1[j] = 0.0;
       v.d2[j] = 0.0;
@@ -270,7 +281,7 @@ TwoPunctures (CCTK_ARGUMENTS)
 
       /* Loop until both ADM masses are within adm_tol of their target */
       do {
-        CCTK_VInfo (CCTK_THORNSTRING, "Bare masses: mp=%g, mm=%g",
+        CCTK_VInfo (CCTK_THORNSTRING, "Bare masses: mp=%.15g, mm=%.15g",
                     (double)*mp, (double)*mm);
         Newton (cctkGH, nvar, n1, n2, n3, v, Newton_tol, 1);
 
@@ -286,8 +297,8 @@ TwoPunctures (CCTK_ARGUMENTS)
         /* 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);
+        CCTK_VInfo (CCTK_THORNSTRING, "ADM mass error: M_p_err=%.15g, M_m_err=%.15g",
+                    (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 */
@@ -299,10 +310,10 @@ TwoPunctures (CCTK_ARGUMENTS)
         
         /* Set the par_m_plus and par_m_minus parameters */
         sprintf (valbuf,"%.17g", (double) *mp);
-        CCTK_ParameterSet ("par_m_plus", "twopunctures", valbuf);
+        CCTK_ParameterSet ("par_m_plus", "TwoPunctures", valbuf);
         
         sprintf (valbuf,"%.17g", (double) *mm);
-        CCTK_ParameterSet ("par_m_minus", "twopunctures", valbuf);
+        CCTK_ParameterSet ("par_m_minus", "TwoPunctures", valbuf);
         
       } while ( (mp_adm_err > adm_tol) ||
                 (mm_adm_err > adm_tol) );
@@ -314,6 +325,8 @@ TwoPunctures (CCTK_ARGUMENTS)
  
     F_of_v (cctkGH, nvar, n1, n2, n3, v, F, u);
 
+    SpecCoef(n1, n2, n3, 0, v.d0, cf_v.d0);
+
     CCTK_VInfo (CCTK_THORNSTRING,
 		  "The two puncture masses are mp=%.17g and mm=%.17g",
                 (double) *mp, (double) *mm);
@@ -456,8 +469,7 @@ TwoPunctures (CCTK_ARGUMENTS)
             (0, nvar, n1, n2, n3, v, xx, yy, zz);
           break;
         case GSM_evaluation:
-          U = PunctIntPolAtArbitPosition
-            (0, nvar, n1, n2, n3, v, xx, yy, zz);
+          U = PunctIntPolAtArbitPositionFast(0, nvar, n1, n2, n3, cf_v, xx, yy, zz);
           break;
         default:
           assert (0);
@@ -664,6 +676,7 @@ TwoPunctures (CCTK_ARGUMENTS)
     free_dvector (F, 0, ntotal - 1);
     free_derivs (&u, ntotal);
     free_derivs (&v, ntotal);
+    free_derivs (&cf_v, ntotal);
   }
 }