Tom's C version of IDAnalyticBH is now the default. Passes testsuites

for Schwarzschild and 2 misner holes, need to test multiple misner
and BL still


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/IDAnalyticBH/trunk@40 6a3ddf76-46e1-4315-99d9-bc56cac1ef84

1 files changed, 219 insertions, 0 deletions

diff --git a/src/Misner_multiple.c b/src/Misner_multiple.c
new file mode 100644
index 0000000..2ba3a1d
--- /dev/null
+++ b/src/Misner_multiple.c
@@ -0,0 +1,219 @@
+ /*@@
+   @file      Misner_multiple.F
+   @date      
+   @author    Carsten Gundlach
+   @desc 
+      Set up initial data for multiple Misner black holes
+   @enddesc 
+ @@*/
+
+#include "cctk.h"
+#include "cctk_arguments.h"
+#include "cctk_parameters.h"
+
+#include "CactusEinstein/Einstein/src/Einstein.h"
+
+#define SQR(a) ((a)*(a))
+
+void Misner_init(int n, CCTK_REAL mu, int terms);
+void MisnerEvalPsi(CCTK_REAL x, CCTK_REAL y, CCTK_REAL z, CCTK_REAL *res); 
+
+ /*@@
+   @routine    Misner_multiple
+   @date       
+   @author     Carsten Gundlach
+   @desc 
+      Set up initial data for multiple Misner black holes
+   @enddesc 
+   @calls     
+   @history 
+ 
+   @endhistory 
+
+@@*/
+void Misner_multiple(CCTK_CARGUMENTS)
+{
+  DECLARE_CCTK_CARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+
+  int i, j, k;
+  CCTK_REAL xval, yval, zval;
+  CCTK_REAL one, zero, tmp0, tmp1, tmp2, tmp3, tmp4;
+  CCTK_REAL nm_eps = 0.000001;   /* finite differencing step*/
+
+  int nx,ny,nz;
+  int index;
+
+  nx = cctk_lsh[0];
+  ny = cctk_lsh[1];
+  nz = cctk_lsh[2];
+
+  one  = 1.0;
+  zero = 0.0;
+
+  /*     Initialize C global variables
+   *     -----------------------------
+   */
+
+  Misner_init(misner_nbh,mu,nmax);
+
+  /*     Get value of psi pointwise from a C function
+   *     --------------------------------------------
+   */
+  for(k=0; k < nz; k++)
+  {
+    for(j=0; j < ny; j++)
+    {
+      for(i=0; i < nx; i++)
+      {
+        index = CCTK_GFINDEX3D(cctkGH, i,j,k);
+
+        xval = x[index];
+        yval = y[index];
+        zval = z[index];
+
+        MisnerEvalPsi(xval,yval,zval,&tmp0);
+        psi[index] = tmp0;
+               
+        /*              Only calculate derivatives of psi if required
+         *              ---------------------------------------------
+         */
+        if (use_conformal_derivs == 1)
+        {
+	  printf("Hello\n");
+          MisnerEvalPsi(xval+nm_eps,yval,zval,&tmp1);
+          MisnerEvalPsi(xval-nm_eps,yval,zval,&tmp2);
+          psix[index] = 0.5*(tmp1-tmp2)/nm_eps;
+          psixx[index] = (tmp1+tmp2-2.0*tmp0)/SQR(nm_eps);
+
+          MisnerEvalPsi(xval,yval+nm_eps,zval,&tmp1);
+          MisnerEvalPsi(xval,yval-nm_eps,zval,&tmp2);
+          psiy[index]  = 0.5*(tmp1-tmp2)/nm_eps;
+          psiyy[index] = (tmp1+tmp2-2.0*tmp0)/SQR(nm_eps);
+
+          MisnerEvalPsi(xval,yval,zval+nm_eps,&tmp1);
+          MisnerEvalPsi(xval,yval,zval-nm_eps,&tmp2);
+          psiz[index] = 0.5*(tmp1-tmp2)/nm_eps;
+          psizz[index] = (tmp1+tmp2-2.0*tmp0)/SQR(nm_eps);
+
+          MisnerEvalPsi(xval+nm_eps,yval+nm_eps,zval,&tmp1);
+          MisnerEvalPsi(xval+nm_eps,yval-nm_eps,zval,&tmp2);
+          MisnerEvalPsi(xval-nm_eps,yval+nm_eps,zval,&tmp3);
+          MisnerEvalPsi(xval-nm_eps,yval-nm_eps,zval,&tmp4);
+          psixy[index] = 0.25*(tmp1-tmp2-tmp3+tmp4)/SQR(nm_eps);
+
+          MisnerEvalPsi(xval,yval+nm_eps,zval+nm_eps,&tmp1);
+          MisnerEvalPsi(xval,yval-nm_eps,zval+nm_eps,&tmp2);
+          MisnerEvalPsi(xval,yval+nm_eps,zval-nm_eps,&tmp3);
+          MisnerEvalPsi(xval,yval-nm_eps,zval-nm_eps,&tmp4);
+          psiyz[index] = 0.25*(tmp1-tmp2-tmp3+tmp4)/SQR(nm_eps);
+
+          MisnerEvalPsi(xval+nm_eps,yval,zval+nm_eps,&tmp1);
+          MisnerEvalPsi(xval+nm_eps,yval,zval-nm_eps,&tmp2);
+          MisnerEvalPsi(xval-nm_eps,yval,zval+nm_eps,&tmp3);
+          MisnerEvalPsi(xval-nm_eps,yval,zval-nm_eps,&tmp4);
+          psixz[index] = 0.25*(tmp1-tmp2-tmp3+tmp4)/SQR(nm_eps);
+
+        }
+
+      }
+    }
+  }
+
+  /*     Cactus conventions
+   *     ------------------
+   */
+
+  if (use_conformal_derivs == 1)
+  {
+    for(k=0; k < nz; k++)
+    {
+      for(j=0; j < ny; j++)
+      {
+        for(i=0; i < nx; i++)
+        {
+          index = CCTK_GFINDEX3D(cctkGH, i,j,k);
+          
+          psix[index]  /= psi[index];
+          psiy[index]  /= psi[index];
+          psiz[index]  /= psi[index];
+          psixx[index] /= psi[index];
+          psixy[index] /= psi[index];
+          psixz[index] /= psi[index];
+          psiyy[index] /= psi[index];
+          psiyz[index] /= psi[index];
+          psizz[index] /= psi[index];
+        }
+      }
+    }
+  }
+
+  /*     Metric depends on conformal state
+   *     ---------------------------------
+   */
+
+  if (*conformal_state == CONFORMAL_METRIC)
+  {
+    for(k=0; k < nz; k++)
+    {
+      for(j=0; j < ny; j++)
+      {
+        for(i=0; i < nx; i++)
+        {
+          index = CCTK_GFINDEX3D(cctkGH, i,j,k);
+          
+          gxx[index] = one;
+          gyy[index] = one;
+          gzz[index] = one;
+          gxy[index] = zero;
+          gxz[index] = zero;
+          gyz[index] = zero;
+        }
+      }
+    }
+  }
+  else 
+  {
+    for(k=0; k < nz; k++)
+    {
+      for(j=0; j < ny; j++)
+      {
+        for(i=0; i < nx; i++)
+        {
+          index = CCTK_GFINDEX3D(cctkGH, i,j,k);
+          
+          gxx[index] = pow(psi[index],4);
+          gyy[index] = gxx[index];
+          gzz[index] = gxx[index];
+          gxy[index] = zero;
+          gxz[index] = zero;
+          gyz[index] = zero;
+        }
+      }
+    }
+  }
+
+  /*     Time-symmetric data
+   *     -------------------
+   */
+
+  for(k=0; k < nz; k++)
+  {
+    for(j=0; j < ny; j++)
+    {
+      for(i=0; i < nx; i++)
+      {
+        index = CCTK_GFINDEX3D(cctkGH, i,j,k);
+          
+        kxx[index] = zero;
+        kyy[index] = zero;
+        kzz[index] = zero;
+        kxy[index] = zero;
+        kxz[index] = zero;
+        kyz[index] = zero;
+      }
+    }
+  }
+
+  return;
+}