implement initial guess, matter rescaling, debug output files

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

3 files changed, 126 insertions, 7 deletions

diff --git a/src/FuncAndJacobian.c b/src/FuncAndJacobian.c
index e105883..87226af 100644
--- a/src/FuncAndJacobian.c
+++ b/src/FuncAndJacobian.c
@@ -218,7 +218,7 @@ F_of_v (CCTK_POINTER_TO_CONST cctkGH,
       for (j = 0; j < n2; j++)
         for (k = 0; k < n3; k++)
         {
-          i3D = i*n2*n3 + j*n3 + k;
+          i3D = Index(0,i,j,k,1,n1,n2,n3);
 
           al = Pih * (2 * i + 1) / n1;
           A = -cos (al);
@@ -260,10 +260,15 @@ F_of_v (CCTK_POINTER_TO_CONST cctkGH,
     for (i = 0; i < n1; i++)
       for (j = 0; j < n2; j++)
         for (k = 0; k < n3; k++)
-          sources[i*n2*n3 + j*n3 + k]=0.0;
+          sources[Index(0,i,j,k,1,n1,n2,n3)]=0.0;
 
   Derivatives_AB3 (nvar, n1, n2, n3, v);
-
+  FILE *debugfile;
+  double psi, psi2, psi4, psi7, r_plus, r_minus;
+  if (do_residuum_debug_output)
+  {
+    debugfile=fopen("res.dat", "w");
+  }
   for (i = 0; i < n1; i++)
   {
     for (j = 0; j < n2; j++)
@@ -301,7 +306,7 @@ F_of_v (CCTK_POINTER_TO_CONST cctkGH,
         // Calculation of (y,z) and
         // (U, U_x, U_y, U_z, U_xx, U_xy, U_xz, U_yy, U_yz, U_zz)
         rx3_To_xyz (nvar, x, r, phi, &y, &z, U);
-        NonLinEquations (sources[i*n2*n3 + j*n3 + k],
+        NonLinEquations (sources[Index(0,i,j,k,1,n1,n2,n3)],
                          A, B, X, R, x, r, phi, y, z, U, values);
         for (ivar = 0; ivar < nvar; ivar++)
         {
@@ -319,8 +324,28 @@ F_of_v (CCTK_POINTER_TO_CONST cctkGH,
           u.d33[indx] = U.d33[ivar];        //      U_zz
         }
       }
+      if (do_residuum_debug_output)
+      {
+        r_plus = sqrt ((x - par_b) * (x - par_b) + y * y + z * z);
+        r_minus = sqrt ((x + par_b) * (x + par_b) + y * y + z * z);
+        psi =
+          1.+ 0.5 * par_m_plus / r_plus + 0.5 * par_m_minus / r_minus + U.d0[0];
+        psi2 = psi * psi;
+        psi4 = psi2 * psi2;
+        psi7 = psi * psi2 * psi4;
+        fprintf(debugfile, "%.8g %.8g %.8g\n", x, y,
+           U.d11[0] +
+           U.d22[0] +
+           U.d33[0]
+           + 2.0 * Pi / psi2/psi * sources[Index(0,i,j,k,1,n1,n2,n3)]
+           );
+      }
     }
   }
+  if (do_residuum_debug_output)
+  {
+    fclose(debugfile);
+  }
   free(sources);
   free_dvector (values, 0, nvar - 1);
   free_derivs (&U, nvar);
diff --git a/src/Newton.c b/src/Newton.c
index 2eefea1..076b926 100644
--- a/src/Newton.c
+++ b/src/Newton.c
@@ -421,9 +421,6 @@ Newton (CCTK_POINTER_TO_CONST cctkGH,
   allocate_derivs (&u, ntotal);
 
 /*         TestRelax(nvar, n1, n2, n3, v, dv.d0); */
-  for (j = 0; j < ntotal; j++)
-    v.d0[j] = 0;
-
   it = 0;
   dmax = 1;
   while (dmax > tol && it < itmax)
@@ -456,6 +453,14 @@ Newton (CCTK_POINTER_TO_CONST cctkGH,
     }
     it += 1;
   }
+  if (itmax==0)
+  {
+      F_of_v (cctkGH, nvar, n1, n2, n3, v, F, u);
+      dmax = -1;
+      for (j = 0; j < ntotal; j++)
+        if (fabs (F[j]) > dmax)
+          dmax = fabs (F[j]);
+  }
   printf ("Newton: it=%d \t |F|=%e \n", it, dmax);
   fflush(stdout);
 
diff --git a/src/TwoPunctures.c b/src/TwoPunctures.c
index 4493e7b..05260ca 100644
--- a/src/TwoPunctures.c
+++ b/src/TwoPunctures.c
@@ -22,6 +22,75 @@ static inline double pow4 (const double x)
   return x*x*x*x;
 }
 
+void set_initial_guess(CCTK_POINTER_TO_CONST cctkGH,
+                       derivs v)
+{
+  DECLARE_CCTK_PARAMETERS;
+
+  int nvar = 1, n1 = npoints_A, n2 = npoints_B, n3 = npoints_phi;
+  CCTK_REAL *s_x, *s_y, *s_z;
+  CCTK_REAL al, A, Am1, be, B, phi, R, r, X;
+  CCTK_INT i, j, k, i3D, ivar, indx;
+  derivs U;
+  FILE *debug_file;
+  CCTK_REAL tmp_r;
+
+  s_x    =calloc(n1*n2*n3, sizeof(CCTK_REAL));
+  s_y    =calloc(n1*n2*n3, sizeof(CCTK_REAL));
+  s_z    =calloc(n1*n2*n3, sizeof(CCTK_REAL));
+  allocate_derivs (&U, nvar);
+  for (i = 0; i < n1; i++)
+    for (j = 0; j < n2; j++)
+      for (k = 0; k < n3; k++)
+      {
+        i3D = Index(0,i,j,k,1,n1,n2,n3);
+
+        al = Pih * (2 * i + 1) / n1;
+        A = -cos (al);
+        be = Pih * (2 * j + 1) / n2;
+        B = -cos (be);
+        phi = 2. * Pi * k / n3;
+
+        // Calculation of (X,R)
+        AB_To_XR (nvar, A, B, &X, &R, U);
+        // Calculation of (x,r)
+        C_To_c (nvar, X, R, &(s_x[i3D]), &r, U);
+        // Calculation of (y,z)
+        rx3_To_xyz (nvar, s_x[i3D], r, phi, &(s_y[i3D]), &(s_z[i3D]), U);
+      }
+  Set_Initial_Guess_for_u(cctkGH, n1*n2*n3, v.d0, s_x, s_y, s_z);
+  
+  if (do_initial_debug_output)
+    debug_file=fopen("initial.dat", "w");
+  for (i = 0; i < n1; i++)
+    for (j = 0; j < n2; j++)
+    {
+      for (k = 0; k < n3; k++)
+      {
+        indx=Index(0,i,j,k,1,n1,n2,n3);
+        tmp_r= s_x[indx]*s_x[indx]+
+               s_y[indx]*s_y[indx]+
+               s_z[indx]*s_z[indx];
+        v.d0[indx]/=-tmp_r/((sqrt(tmp_r)+10)*(sqrt(tmp_r)+10))+1;
+      }
+      indx=Index(0,i,j,0,1,n1,n2,n3);
+      if (do_initial_debug_output)
+      {
+        fprintf(debug_file, "%.8g %.8g %.8g %.8g %.8g\n", s_x[indx], s_y[indx],
+                v.d0[indx],
+                v.d0[indx]*(-cos(Pih * (2 * i + 1) / n1)-1.0),
+                (-cos(Pih * (2 * i + 1) / n1)-1.0));
+      }
+    }
+  if (do_initial_debug_output)
+    fclose(debug_file);
+  Derivatives_AB3 (nvar, n1, n2, n3, v);
+  free(s_z);
+  free(s_y);
+  free(s_x);
+  free_derivs (&U, nvar);
+}
+
 // -------------------------------------------------------------------
 void
 TwoPunctures (CCTK_ARGUMENTS)
@@ -48,6 +117,16 @@ TwoPunctures (CCTK_ARGUMENTS)
     allocate_derivs (&v, ntotal);
 
     CCTK_INFO ("Solving puncture equation");
+    /* initialise to 0 */
+    for (j = 0; j < ntotal; j++)
+      v.d0[j] = 0.0;
+    /* call for external initial guess */
+    if (use_external_initial_guess)
+    {
+      set_initial_guess(cctkGH, v);
+      F_of_v (cctkGH, nvar, n1, n2, n3, v, F, u);
+    }
+
     Newton (cctkGH, nvar, n1, n2, n3, v, Newton_tol, Newton_maxit);
 
     F_of_v (cctkGH, nvar, n1, n2, n3, v, F, u);
@@ -237,6 +316,16 @@ TwoPunctures (CCTK_ARGUMENTS)
     }
   }
 
+  if (use_sources && rescale_sources)
+  {
+    Rescale_Sources(cctkGH,
+                    cctk_lsh[0]*cctk_lsh[1]*cctk_lsh[2],
+                    x, y, z,
+                    psi,
+                    gxx, gyy, gzz,
+                    gxy, gxz, gyz);
+  }
+
   if (0) {
     /* Keep the result around for the next time */
     free_dvector (F, 0, ntotal - 1);