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

1 files changed, 29 insertions, 4 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);