- add source terms provided by aliased function

- small fix suggested by Marcus concerning the initialisation of the dvs


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

10 files changed, 201 insertions, 79 deletions

diff --git a/interface.ccl b/interface.ccl
index a969173..3bb592b 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -5,6 +5,14 @@ IMPLEMENTS: TwoPunctures
 
 INHERITS: ADMBase StaticConformal grid
 
-
-
 REAL puncture_u TYPE=gf
+
+CCTK_INT FUNCTION Rho_ADM(                 \
+    CCTK_POINTER_TO_CONST IN cctkGH,       \
+    CCTK_INT              IN size,         \
+    CCTK_POINTER          IN source,       \
+    CCTK_POINTER          IN x,            \
+    CCTK_POINTER          IN y,            \
+    CCTK_POINTER          IN z             \
+    )
+USES FUNCTION Rho_ADM
diff --git a/param.ccl b/param.ccl
index 3703c97..82463f6 100644
--- a/param.ccl
+++ b/param.ccl
@@ -101,3 +101,10 @@ REAL par_S_minus[3] "spin of the m- puncture"
 {
   (*:*) :: ""
 } 0.0
+
+
+INT par_use_sources "Use sources?"
+{
+  0:1 :: "0 for no (default), 1 for yes"
+} 0
+
diff --git a/schedule.ccl b/schedule.ccl
index 797d4d6..37262fb 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -5,7 +5,10 @@ if (keep_u_around) {
   STORAGE: puncture_u
 }
 
-
+SCHEDULE TwoPunctures_ParamCheck AT PARAMCHECK
+{
+  LANG: C
+} "Check parameters and thorn needs"
 
 SCHEDULE TwoPunctures IN ADMBase_InitialData
 {
diff --git a/src/Equations.c b/src/Equations.c
index e1de7b5..0ea280a 100644
--- a/src/Equations.c
+++ b/src/Equations.c
@@ -136,7 +136,8 @@ BY_Aijofxyz (double x, double y, double z, double Aij[3][3])
 //*******           Nonlinear Equations                **********
 //---------------------------------------------------------------
 void
-NonLinEquations (double A, double B, double X, double R,
+NonLinEquations (CCTK_REAL rho_adm,
+     double A, double B, double X, double R,
 		 double x, double r, double phi,
 		 double y, double z, derivs U, double *values)
 {
@@ -154,7 +155,8 @@ NonLinEquations (double A, double B, double X, double R,
   psi7 = psi * psi2 * psi4;
 
   values[0] =
-    U.d11[0] + U.d22[0] + U.d33[0] + 0.125 * BY_KKofxyz (x, y, z) / psi7;
+    U.d11[0] + U.d22[0] + U.d33[0] + 0.125 * BY_KKofxyz (x, y, z) / psi7 +
+    2.0 * Pi * psi4*psi * rho_adm;
 
 }
 
diff --git a/src/FuncAndJacobian.c b/src/FuncAndJacobian.c
index 7c30da0..c2c972e 100644
--- a/src/FuncAndJacobian.c
+++ b/src/FuncAndJacobian.c
@@ -7,6 +7,7 @@
 #include <math.h>
 #include <ctype.h>
 #include <time.h>
+#include "cctk.h"
 #include "cctk_Parameters.h"
 #include "TP_utilities.h"
 #include "TwoPunctures.h"
@@ -188,19 +189,78 @@ Derivatives_AB3 (int nvar, int n1, int n2, int n3, derivs v)
 
 // -------------------------------------------------------------------------------
 void
-F_of_v (int nvar, int n1, int n2, int n3, derivs v, double *F,
-	derivs u)
+F_of_v (CCTK_POINTER_TO_CONST cctkGH,
+        int nvar, int n1, int n2, int n3, derivs v, double *F,
+        derivs u)
 {				//      Calculates the left hand sides of the non-linear equations F_m(v_n)=0
   //      and the function u (u.d0[]) as well as its derivatives
   //      (u.d1[], u.d2[], u.d3[], u.d11[], u.d12[], u.d13[], u.d22[], u.d23[], u.d33[])
   //      at interior points and at the boundaries "+/-"
+  DECLARE_CCTK_PARAMETERS;
   int i, j, k, ivar, indx;
   double al, be, A, B, X, R, x, r, phi, y, z, Am1, *values;
   derivs U;
+  CCTK_REAL *sources;
 
   values = dvector (0, nvar - 1);
   allocate_derivs (&U, nvar);
 
+  if (par_use_sources)
+  {
+    CCTK_REAL *s_x, *s_y, *s_z;
+    CCTK_INT i3D;
+    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));
+    sources=calloc(n1*n2*n3, sizeof(CCTK_REAL));
+    for (i = 0; i < n1; i++)
+      for (j = 0; j < n2; j++)
+        for (k = 0; k < n3; k++)
+        {
+          i3D = i*n2*n3 + j*n3 + k;
+
+          al = Pih * (2 * i + 1) / n1;
+          A = -cos (al);
+          be = Pih * (2 * j + 1) / n2;
+          B = -cos (be);
+          phi = 2. * Pi * k / n3;
+
+          Am1 = A - 1;
+          for (ivar = 0; ivar < nvar; ivar++)
+          {
+            indx = Index (ivar, i, j, k, nvar, n1, n2, n3);
+            U.d0[ivar] = Am1 * v.d0[indx];        // U
+            U.d1[ivar] = v.d0[indx] + Am1 * v.d1[indx];        // U_A
+            U.d2[ivar] = Am1 * v.d2[indx];        // U_B
+            U.d3[ivar] = Am1 * v.d3[indx];        // U_3
+            U.d11[ivar] = 2 * v.d1[indx] + Am1 * v.d11[indx];        // U_AA
+            U.d12[ivar] = v.d2[indx] + Am1 * v.d12[indx];        // U_AB
+            U.d13[ivar] = v.d3[indx] + Am1 * v.d13[indx];        // U_AB
+            U.d22[ivar] = Am1 * v.d22[indx];        // U_BB
+            U.d23[ivar] = Am1 * v.d23[indx];        // U_B3
+            U.d33[ivar] = Am1 * v.d33[indx];        // U_33
+          }
+          // Calculation of (X,R) and
+          // (U_X, U_R, U_3, U_XX, U_XR, U_X3, U_RR, U_R3, U_33)
+          AB_To_XR (nvar, A, B, &X, &R, U);
+          // Calculation of (x,r) and
+          // (U, U_x, U_r, U_3, U_xx, U_xr, U_x3, U_rr, U_r3, U_33)
+          C_To_c (nvar, X, R, &(s_x[i3D]), &r, U);
+          // 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, s_x[i3D], r, phi, &(s_y[i3D]), &(s_z[i3D]), U);
+        }
+    Rho_ADM(cctkGH, n1*n2*n3, sources, s_x, s_y, s_z);
+    free(s_z);
+    free(s_y);
+    free(s_x);
+  }
+  else
+    for (i = 0; i < n1; i++)
+      for (j = 0; j < n2; j++)
+        for (k = 0; k < n3; k++)
+          sources[i*n1*n2 + j*n2 + k]=0.0;
+
   Derivatives_AB3 (nvar, n1, n2, n3, v);
 
   for (i = 0; i < n1; i++)
@@ -210,55 +270,58 @@ F_of_v (int nvar, int n1, int n2, int n3, derivs v, double *F,
       for (k = 0; k < n3; k++)
       {
 
-	al = Pih * (2 * i + 1) / n1;
-	A = -cos (al);
-	be = Pih * (2 * j + 1) / n2;
-	B = -cos (be);
-	phi = 2. * Pi * k / n3;
-
-	Am1 = A - 1;
-	for (ivar = 0; ivar < nvar; ivar++)
-	{
-	  indx = Index (ivar, i, j, k, nvar, n1, n2, n3);
-	  U.d0[ivar] = Am1 * v.d0[indx];	// U
-	  U.d1[ivar] = v.d0[indx] + Am1 * v.d1[indx];	// U_A
-	  U.d2[ivar] = Am1 * v.d2[indx];	// U_B
-	  U.d3[ivar] = Am1 * v.d3[indx];	// U_3
-	  U.d11[ivar] = 2 * v.d1[indx] + Am1 * v.d11[indx];	// U_AA
-	  U.d12[ivar] = v.d2[indx] + Am1 * v.d12[indx];	// U_AB
-	  U.d13[ivar] = v.d3[indx] + Am1 * v.d13[indx];	// U_AB
-	  U.d22[ivar] = Am1 * v.d22[indx];	// U_BB
-	  U.d23[ivar] = Am1 * v.d23[indx];	// U_B3
-	  U.d33[ivar] = Am1 * v.d33[indx];	// U_33
-	}
-	// Calculation of (X,R) and
-	// (U_X, U_R, U_3, U_XX, U_XR, U_X3, U_RR, U_R3, U_33)
-	AB_To_XR (nvar, A, B, &X, &R, U);
-	// Calculation of (x,r) and
-	// (U, U_x, U_r, U_3, U_xx, U_xr, U_x3, U_rr, U_r3, U_33)
-	C_To_c (nvar, X, R, &x, &r, U);
-	// 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 (A, B, X, R, x, r, phi, y, z, U, values);
-	for (ivar = 0; ivar < nvar; ivar++)
-	{
-	  indx = Index (ivar, i, j, k, nvar, n1, n2, n3);
-	  F[indx] = values[ivar] * FAC;
-	  u.d0[indx] = U.d0[ivar];	//  U
-	  u.d1[indx] = U.d1[ivar];	//      U_x
-	  u.d2[indx] = U.d2[ivar];	//      U_y
-	  u.d3[indx] = U.d3[ivar];	//      U_z
-	  u.d11[indx] = U.d11[ivar];	//      U_xx
-	  u.d12[indx] = U.d12[ivar];	//      U_xy
-	  u.d13[indx] = U.d13[ivar];	//      U_xz
-	  u.d22[indx] = U.d22[ivar];	//      U_yy
-	  u.d23[indx] = U.d23[ivar];	//      U_yz
-	  u.d33[indx] = U.d33[ivar];	//      U_zz
-	}
+        al = Pih * (2 * i + 1) / n1;
+        A = -cos (al);
+        be = Pih * (2 * j + 1) / n2;
+        B = -cos (be);
+        phi = 2. * Pi * k / n3;
+
+        Am1 = A - 1;
+        for (ivar = 0; ivar < nvar; ivar++)
+        {
+          indx = Index (ivar, i, j, k, nvar, n1, n2, n3);
+          U.d0[ivar] = Am1 * v.d0[indx];        // U
+          U.d1[ivar] = v.d0[indx] + Am1 * v.d1[indx];        // U_A
+          U.d2[ivar] = Am1 * v.d2[indx];        // U_B
+          U.d3[ivar] = Am1 * v.d3[indx];        // U_3
+          U.d11[ivar] = 2 * v.d1[indx] + Am1 * v.d11[indx];        // U_AA
+          U.d12[ivar] = v.d2[indx] + Am1 * v.d12[indx];        // U_AB
+          U.d13[ivar] = v.d3[indx] + Am1 * v.d13[indx];        // U_AB
+          U.d22[ivar] = Am1 * v.d22[indx];        // U_BB
+          U.d23[ivar] = Am1 * v.d23[indx];        // U_B3
+          U.d33[ivar] = Am1 * v.d33[indx];        // U_33
+        }
+        // Calculation of (X,R) and
+        // (U_X, U_R, U_3, U_XX, U_XR, U_X3, U_RR, U_R3, U_33)
+        AB_To_XR (nvar, A, B, &X, &R, U);
+        // Calculation of (x,r) and
+        // (U, U_x, U_r, U_3, U_xx, U_xr, U_x3, U_rr, U_r3, U_33)
+        C_To_c (nvar, X, R, &x, &r, U);
+        // 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],
+                         A, B, X, R, x, r, phi, y, z, U, values);
+        for (ivar = 0; ivar < nvar; ivar++)
+        {
+          indx = Index (ivar, i, j, k, nvar, n1, n2, n3);
+          F[indx] = values[ivar] * FAC;
+          u.d0[indx] = U.d0[ivar];        //  U
+          u.d1[indx] = U.d1[ivar];        //      U_x
+          u.d2[indx] = U.d2[ivar];        //      U_y
+          u.d3[indx] = U.d3[ivar];        //      U_z
+          u.d11[indx] = U.d11[ivar];        //      U_xx
+          u.d12[indx] = U.d12[ivar];        //      U_xy
+          u.d13[indx] = U.d13[ivar];        //      U_xz
+          u.d22[indx] = U.d22[ivar];        //      U_yy
+          u.d23[indx] = U.d23[ivar];        //      U_yz
+          u.d33[indx] = U.d33[ivar];        //      U_zz
+        }
       }
     }
   }
+  if (par_use_sources)
+    free(sources);
   free_dvector (values, 0, nvar - 1);
   free_derivs (&U, nvar);
 }
diff --git a/src/Newton.c b/src/Newton.c
index 593789f..0dbb627 100644
--- a/src/Newton.c
+++ b/src/Newton.c
@@ -173,7 +173,8 @@ relax (double *dv, int nvar, int n1, int n2, int n3,
 
 // -----------------------------------------------------------------------------------
 void
-TestRelax (int nvar, int n1, int n2, int n3, derivs v,
+TestRelax (CCTK_POINTER_TO_CONST cctkGH,
+     int nvar, int n1, int n2, int n3, derivs v,
 	   double *dv)
 {
   DECLARE_CCTK_PARAMETERS;
@@ -190,7 +191,7 @@ TestRelax (int nvar, int n1, int n2, int n3, derivs v,
   cols = imatrix (0, ntotal - 1, 0, maxcol - 1);
   ncols = ivector (0, ntotal - 1);
 
-  F_of_v (nvar, n1, n2, n3, v, F, u);
+  F_of_v (cctkGH, nvar, n1, n2, n3, v, F, u);
 
   SetMatrix_JFD (nvar, n1, n2, n3, u, ncols, cols, JFD);
 
@@ -225,8 +226,9 @@ TestRelax (int nvar, int n1, int n2, int n3, derivs v,
 
 // -----------------------------------------------------------------------------------
 int
-bicgstab (int nvar, int n1, int n2, int n3, derivs v,
-	  derivs dv, int output, int itmax, double tol, double *normres)
+bicgstab (CCTK_POINTER_TO_CONST cctkGH,
+          int nvar, int n1, int n2, int n3, derivs v,
+	        derivs dv, int output, int itmax, double tol, double *normres)
 {
   DECLARE_CCTK_PARAMETERS;
   int ntotal = n1 * n2 * n3 * nvar, ii, j;
@@ -246,7 +248,7 @@ bicgstab (int nvar, int n1, int n2, int n3, derivs v,
   cols = imatrix (0, ntotal - 1, 0, maxcol - 1);
   ncols = ivector (0, ntotal - 1);
 
-  F_of_v (nvar, n1, n2, n3, v, F, u);
+  F_of_v (cctkGH, nvar, n1, n2, n3, v, F, u);
   SetMatrix_JFD (nvar, n1, n2, n3, u, ncols, cols, JFD);
 
   /* temporary storage */
@@ -389,7 +391,8 @@ bicgstab (int nvar, int n1, int n2, int n3, derivs v,
 
 // -------------------------------------------------------------------
 void
-Newton (int nvar, int n1, int n2, int n3,
+Newton (CCTK_POINTER_TO_CONST cctkGH,
+  int nvar, int n1, int n2, int n3,
 	derivs v, double tol, int itmax)
 {
   DECLARE_CCTK_PARAMETERS;
@@ -401,7 +404,7 @@ Newton (int nvar, int n1, int n2, int n3,
   allocate_derivs (&dv, ntotal);
   allocate_derivs (&u, ntotal);
 
-/* 	TestRelax(nvar, n1, n2, n3, v, dv.d0); */
+/*         TestRelax(nvar, n1, n2, n3, v, dv.d0); */
   for (j = 0; j < ntotal; j++)
     v.d0[j] = 0;
 
@@ -411,28 +414,28 @@ Newton (int nvar, int n1, int n2, int n3,
   {
     if (it == 0)
     {
-      F_of_v (nvar, n1, n2, n3, v, F, u);
+      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]);
-	dv.d0[j] = 0;
-      }
+        if (fabs (F[j]) > dmax)
+          dmax = fabs (F[j]);
     }
+    for (j = 0; j < ntotal; j++)
+      dv.d0[j] = 0;
     printf ("Newton: it=%d \t |F|=%e\n", it, dmax);
     fflush(stdout);
     ii =
-      bicgstab (nvar, n1, n2, n3, v, dv, verbose, 100, dmax * 1.e-3, &normres);
+      bicgstab (cctkGH,
+                nvar, n1, n2, n3, v, dv, verbose, 100, dmax * 1.e-3, &normres);
     for (j = 0; j < ntotal; j++)
       v.d0[j] -= dv.d0[j];
-    F_of_v (nvar, n1, n2, n3, v, F, u);
+    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]);
+        dmax = fabs (F[j]);
       }
     }
     it += 1;
diff --git a/src/ParamCheck.c b/src/ParamCheck.c
new file mode 100644
index 0000000..ca06f64
--- /dev/null
+++ b/src/ParamCheck.c
@@ -0,0 +1,32 @@
+// TwoPunctures:  File  "TwoPunctures.c"
+
+#include <assert.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <ctype.h>
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+#include "TP_utilities.h"
+#include "TwoPunctures.h"
+
+// -------------------------------------------------------------------
+void
+TwoPunctures_ParamCheck (CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  if (par_use_sources)
+  {
+    CCTK_INFO("Solving for BH-NS");
+    if (CCTK_IsFunctionAliased ("Rho_ADM"))
+      CCTK_INFO("Aliased Functions found");
+    else
+      CCTK_WARN(0, "I found no (aliased) function for matter sources, but "
+                   "was said to use matter.\n");
+  }
+  else
+    CCTK_INFO("not using sources (only BHs)");
+}
diff --git a/src/TwoPunctures.c b/src/TwoPunctures.c
index 1c792b6..e926b65 100644
--- a/src/TwoPunctures.c
+++ b/src/TwoPunctures.c
@@ -46,9 +46,9 @@ TwoPunctures (CCTK_ARGUMENTS)
     allocate_derivs (&v, ntotal);
 
     CCTK_INFO ("Solving puncture equation");
-    Newton (nvar, n1, n2, n3, v, Newton_tol, Newton_maxit);
+    Newton (cctkGH, nvar, n1, n2, n3, v, Newton_tol, Newton_maxit);
 
-    F_of_v (nvar, n1, n2, n3, v, F, u);
+    F_of_v (cctkGH, nvar, n1, n2, n3, v, F, u);
     
     /* print out ADM mass, eq.: \Delta M_ADM=2*r*u=4*b*V for A=1,B=0,phi=0 */
     admMass = (par_m_plus + par_m_minus
diff --git a/src/TwoPunctures.h b/src/TwoPunctures.h
index 9915181..19763e7 100644
--- a/src/TwoPunctures.h
+++ b/src/TwoPunctures.h
@@ -30,7 +30,8 @@ int Index (int ivar, int i, int j, int k, int nvar, int n1, int n2, int n3);
 void allocate_derivs (derivs * v, int n);
 void free_derivs (derivs * v, int n);
 void Derivatives_AB3 (int nvar, int n1, int n2, int n3, derivs v);
-void F_of_v (int nvar, int n1, int n2, int n3, derivs v,
+void F_of_v (CCTK_POINTER_TO_CONST cctkGH,
+       int nvar, int n1, int n2, int n3, derivs v,
 	     double *F, derivs u);
 void J_times_dv (int nvar, int n1, int n2, int n3, derivs dv,
 		 double *Jdv, derivs u);
@@ -61,7 +62,8 @@ void rx3_To_xyz (int nvar, double x, double r, double phi, double *y,
 // Routines in  "Equations.c"
 double BY_KKofxyz (double x, double y, double z);
 void BY_Aijofxyz (double x, double y, double z, double Aij[3][3]);
-void NonLinEquations (double A, double B, double X, double R,
+void NonLinEquations (CCTK_REAL rho_adm,
+          double A, double B, double X, double R,
 		      double x, double r, double phi,
 		      double y, double z, derivs U, double *values);
 void LinEquations (double A, double B, double X, double R,
@@ -77,12 +79,14 @@ void LineRelax_be (double *dv, int i, int k, int nvar, int n1, int n2, int n3,
 		   double *rhs, int *ncols, int **cols, double **JFD);
 void relax (double *dv, int nvar, int n1, int n2, int n3, double *rhs,
 	    int *ncols, int **cols, double **JFD);
-void TestRelax (int nvar, int n1, int n2, int n3, derivs v,
-		double *dv);
-int bicgstab (int nvar, int n1, int n2, int n3, derivs v,
+void TestRelax (CCTK_POINTER_TO_CONST cctkGH,
+                int nvar, int n1, int n2, int n3, derivs v, double *dv);
+int bicgstab (CCTK_POINTER_TO_CONST cctkGH,
+              int nvar, int n1, int n2, int n3, derivs v,
 	      derivs dv, int output, int itmax, double tol, double *normres);
-void Newton (int nvar, int n1, int n2, int n3, derivs v,
-	     double tol, int itmax);
+void Newton (CCTK_POINTER_TO_CONST cctkGH,
+             int nvar, int n1, int n2, int n3, derivs v,
+	           double tol, int itmax);
 
 
 /* 
diff --git a/src/make.code.defn b/src/make.code.defn
index 726aec8..512024e 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -2,7 +2,7 @@
 # $Header$
 
 # Source files in this directory
-SRCS = CoordTransf.c Equations.c FuncAndJacobian.c Newton.c TwoPunctures.c TP_utilities.c
+SRCS = CoordTransf.c Equations.c FuncAndJacobian.c Newton.c TwoPunctures.c TP_utilities.c ParamCheck.c
 
 # Subdirectories containing source files
 SUBDIRS =