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, 143 insertions, 0 deletions

diff --git a/src/FuncAndJacobian.c b/src/FuncAndJacobian.c
index 294254d..8703aa1 100644
--- a/src/FuncAndJacobian.c
+++ b/src/FuncAndJacobian.c
@@ -841,4 +841,147 @@ PunctIntPolAtArbitPosition (int ivar, int nvar, int n1,
   return Ui;
 }
 
+
+//////////////////////////////////////////////////////
+/// Fast Spectral Interpolation Routine Stuff
+//////////////////////////////////////////////////////
+
+
+/* Calculates the value of v at an arbitrary position (A,B,phi)* using the fast routine */
+CCTK_REAL 
+PunctEvalAtArbitPositionFast (CCTK_REAL *v, int ivar, CCTK_REAL A, CCTK_REAL B, CCTK_REAL phi, int nvar, int n1, int n2, int n3)
+{
+  int i, j, k, N;
+  CCTK_REAL *p, *values1, **values2, result;
+  // VASILIS: Nothing should be changed in this routine. This is used by PunctIntPolAtArbitPositionFast
+
+  N = maximum3 (n1, n2, n3);
+  
+  p = dvector (0, N);
+  values1 = dvector (0, N);
+  values2 = dmatrix (0, N, 0, N);
+
+  for (k = 0; k < n3; k++)
+  {
+    for (j = 0; j < n2; j++)
+    {
+      for (i = 0; i < n1; i++) p[i] = v[ivar + nvar * (i + n1 * (j + n2 * k))];
+      //      chebft_Zeros (p, n1, 0);
+      values2[j][k] = chebev (-1, 1, p, n1, A);
+    }
+  }
+
+  for (k = 0; k < n3; k++)
+  {
+    for (j = 0; j < n2; j++) p[j] = values2[j][k];
+    //    chebft_Zeros (p, n2, 0);
+    values1[k] = chebev (-1, 1, p, n2, B);
+  }
+
+  //  fourft (values1, n3, 0);
+  result = fourev (values1, n3, phi);
+  
+  free_dvector (p, 0, N);
+  free_dvector (values1, 0, N);
+  free_dmatrix (values2, 0, N, 0, N);
+
+  return result;
+  //  */
+  //  return 0.;
+}
+
+
+// --------------------------------------------------------------------------*/
+// Calculates the value of v at an arbitrary position (x,y,z) if the spectral coefficients are known //
 /* --------------------------------------------------------------------------*/
+CCTK_REAL
+PunctIntPolAtArbitPositionFast (int ivar, int nvar, int n1,
+			    int n2, int n3, derivs v, CCTK_REAL x, CCTK_REAL y,
+			    CCTK_REAL z)
+{
+  DECLARE_CCTK_PARAMETERS;
+  CCTK_REAL xs, ys, zs, rs2, phi, X, R, A, B, aux1, aux2, result, Ui;
+  // VASILIS: Here the struct derivs v refers to the spectral coeffiecients of variable v not the variable v itself
+
+  xs = x / par_b;
+  ys = y / par_b;
+  zs = z / par_b;
+  rs2 = ys * ys + zs * zs;
+  phi = atan2 (z, y);
+  if (phi < 0)
+    phi += 2 * Pi;
+
+  aux1 = 0.5 * (xs * xs + rs2 - 1);
+  aux2 = sqrt (aux1 * aux1 + rs2);
+  X = asinh (sqrt (aux1 + aux2));
+  R = asin (min(1.0, sqrt (-aux1 + aux2)));
+  if (x < 0)
+    R = Pi - R;
+
+  A = 2 * tanh (0.5 * X) - 1;
+  B = tan (0.5 * R - Piq);
+
+  result = PunctEvalAtArbitPositionFast (v.d0, ivar, A, B, phi, nvar, n1, n2, n3);
+
+  Ui = (A - 1) * result;
+
+  return Ui;
+}
+
+// Evaluates the spectral expansion coefficients of v  
+void SpecCoef(int n1, int n2, int n3, int ivar, CCTK_REAL *v, CCTK_REAL *cf)
+{
+  DECLARE_CCTK_PARAMETERS;
+  // VASILIS: Here v is a pointer to the values of the variable v at the collocation points and cf_v a pointer to the spectral coefficients that this routine calculates
+
+	int i, j, k, N, n, l, m;
+	double *p, ***values3, ***values4;
+	
+	N=maximum3(n1,n2,n3);
+	p=dvector(0,N);
+	values3=d3tensor(0,n1,0,n2,0,n3);
+	values4=d3tensor(0,n1,0,n2,0,n3);
+
+
+
+	      // Caclulate values3[n,j,k] = a_n^{j,k} = (sum_i^(n1-1) f(A_i,B_j,phi_k) Tn(-A_i))/k_n , k_n = N/2 or N 
+	      for(k=0;k<n3;k++) {
+		for(j=0;j<n2;j++) {
+
+		  for(i=0;i<n1;i++) p[i]=v[ivar + (i + n1 * (j + n2 * k))];
+		  
+		  chebft_Zeros(p,n1,0); 
+		  for (n=0;n<n1;n++)	{
+		    values3[n][j][k] = p[n]; 
+		  }
+		}
+	      }
+	    
+	      // Caclulate values4[n,l,k] = a_{n,l}^{k} = (sum_j^(n2-1) a_n^{j,k} Tn(B_j))/k_l , k_l = N/2 or N 
+
+	      for (n = 0; n < n1; n++){
+		for(k=0;k<n3;k++) {
+		  for(j=0;j<n2;j++) p[j]=values3[n][j][k];
+		  chebft_Zeros(p,n2,0);   
+		  for (l = 0; l < n2; l++){
+		  values4[n][l][k] = p[l];
+		  }
+		}
+	      }
+
+	      // Caclulate coefficients  a_{n,l,m} = (sum_k^(n3-1) a_{n,m}^{k} fourier(phi_k))/k_m , k_m = N/2 or N 
+	      for (i = 0; i < n1; i++){
+		for (j = 0; j < n2; j++){
+		  for(k=0;k<n3;k++) p[k]=values4[i][j][k];
+		  fourft(p,n3,0); 
+		  for (k = 0; k<n3; k++){
+		    cf[ivar + (i + n1 * (j + n2 * k))] = p[k];
+		  }
+		}
+	      }
+   
+	      free_dvector(p,0,N);
+	      free_d3tensor(values3,0,n1,0,n2,0,n3);
+	      free_d3tensor(values4,0,n1,0,n2,0,n3);
+
+}