Some fixes, some cleanup.

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

8 files changed, 226 insertions, 75 deletions

diff --git a/param.ccl b/param.ccl
index e16e785..ca2fc22 100644
--- a/param.ccl
+++ b/param.ccl
@@ -47,6 +47,18 @@ INT npoints_phi "Number of coefficients in the phi direction"
 
 
 
+REAL Newton_tol "Tolerance for Newton solver"
+{
+  (0:*) :: ""
+} 1.0e-10
+
+INT Newton_maxit "Maximum number of Newton iterations"
+{
+  0:* :: ""
+} 5
+
+
+
 REAL par_b "x coordinate of the m+ puncture"
 {
   (0.0:*) :: ""
diff --git a/src/Equations.c b/src/Equations.c
index b0d426e..e1de7b5 100644
--- a/src/Equations.c
+++ b/src/Equations.c
@@ -67,8 +67,8 @@ BY_KKofxyz (double x, double y, double z)
                  + np_Pp * n_plus[i] * n_plus[j]) / r2_plus
 	+ 1.5 * (par_P_minus[i] * n_minus[j] + par_P_minus[j] * n_minus[i]
 		 + nm_Pm * n_minus[i] * n_minus[j]) / r2_minus
-	- 3.0 * (np_Sp[i] * n_plus[j] + np_Sp[j] * n_plus[j]) / r3_plus
-	- 3.0 * (nm_Sm[i] * n_minus[j] + nm_Sm[i] * n_minus[j]) / r3_minus;
+	- 3.0 * (np_Sp[i] * n_plus[j] + np_Sp[j] * n_plus[i]) / r3_plus
+	- 3.0 * (nm_Sm[i] * n_minus[j] + nm_Sm[j] * n_minus[i]) / r3_minus;
       if (i == j)
 	Aij -= +1.5 * (np_Pp / r2_plus + nm_Pm / r2_minus);
       AijAij += Aij * Aij;
@@ -124,8 +124,8 @@ BY_Aijofxyz (double x, double y, double z, double Aij[3][3])
 		 + np_Pp * n_plus[i] * n_plus[j]) / r2_plus
 	+ 1.5 * (par_P_minus[i] * n_minus[j] + par_P_minus[j] * n_minus[i]
 		 + nm_Pm * n_minus[i] * n_minus[j]) / r2_minus
-	- 3.0 * (np_Sp[i] * n_plus[j] + np_Sp[j] * n_plus[j]) / r3_plus
-	- 3.0 * (nm_Sm[i] * n_minus[j] + nm_Sm[i] * n_minus[j]) / r3_minus;
+	- 3.0 * (np_Sp[i] * n_plus[j] + np_Sp[j] * n_plus[i]) / r3_plus
+	- 3.0 * (nm_Sm[i] * n_minus[j] + nm_Sm[j] * n_minus[i]) / r3_minus;
       if (i == j)
 	Aij[i][j] -= +1.5 * (np_Pp / r2_plus + nm_Pm / r2_minus);
     }
diff --git a/src/FuncAndJacobian.c b/src/FuncAndJacobian.c
index 568ac13..96adeb4 100644
--- a/src/FuncAndJacobian.c
+++ b/src/FuncAndJacobian.c
@@ -1,5 +1,6 @@
 // TwoPunctures:  File  "FuncAndJacobian.c"
 
+#include <assert.h>
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
@@ -14,6 +15,11 @@
 //#define FAC sin(al)*sin(be)*sin(al)*sin(be)
 //#define FAC 1
 
+static inline double min (double const x, double const y)
+{
+  return x<y ? x : y;
+}
+
 // -------------------------------------------------------------------------------
 int
 Index (int ivar, int i, int j, int k, int nvar, int n1, int n2, int n3)
@@ -648,15 +654,15 @@ interpol (double a, double b, double c, derivs v)
 }
 
 // -------------------------------------------------------------------------------
-// Calculates the value of v at an arbitrary position (A,B,phi)
+// Calculates the value of v at an arbitrary position (x,y,z)
 double
-PunctIntPolAtArbitPosition (int ivar, int nvar, int n1,
-			    int n2, int n3, derivs v, double x, double y,
-			    double z)
+PunctTaylorExpandAtArbitPosition (int ivar, int nvar, int n1,
+                                  int n2, int n3, derivs v, double x, double y,
+                                  double z)
 {
   DECLARE_CCTK_PARAMETERS;
   double xs, ys, zs, rs2, phi, X, R, A, B, al, be, aux1, aux2, a, b, c,
-    result, Us, Ui;
+    result, Ui;
   int i, j, k;
   derivs vv;
   allocate_derivs (&vv, 1);
@@ -672,7 +678,7 @@ PunctIntPolAtArbitPosition (int ivar, int nvar, int n1,
   aux1 = 0.5 * (xs * xs + rs2 - 1);
   aux2 = sqrt (aux1 * aux1 + rs2);
   X = asinh (sqrt (aux1 + aux2));
-  R = asin (sqrt (-aux1 + aux2));
+  R = asin (min(1.0, sqrt (-aux1 + aux2)));
   if (x < 0)
     R = Pi - R;
 
@@ -684,7 +690,7 @@ PunctIntPolAtArbitPosition (int ivar, int nvar, int n1,
   i = rint (al * n1 / Pi - 0.5);
   j = rint (be * n2 / Pi - 0.5);
   k = rint (0.5 * phi * n3 / Pi);
-
+  
   a = al - Pi * (i + 0.5) / n1;
   b = be - Pi * (j + 0.5) / n2;
   c = phi - 2 * Pi * k / n3;
@@ -693,52 +699,44 @@ PunctIntPolAtArbitPosition (int ivar, int nvar, int n1,
   result = interpol (a, b, c, vv);
   free_derivs (&vv, 1);
 
-//      Us = (A-1)*PunctEvalAtArbitPosition(v.d0, A, B, phi, n1, n2, n3);
   Ui = (A - 1) * result;
-//      printf("%e %e %e Us=%e Ui=%e 1-Ui/Us=%e\n",x,y,z,Us,Ui,1-Ui/Us);
 
   return Ui;
+}
 
-/*	calculate_derivs(i+1, j,   k,   ivar, nvar, n1, n2, n3, v, vv, 2);
-	calculate_derivs(i,   j+1, k,   ivar, nvar, n1, n2, n3, v, vv, 3);
-	calculate_derivs(i+1, j+1, k,   ivar, nvar, n1, n2, n3, v, vv, 4);
-	calculate_derivs(i,   j,   k+1, ivar, nvar, n1, n2, n3, v, vv, 5);
-	calculate_derivs(i+1, j,   k+1, ivar, nvar, n1, n2, n3, v, vv, 6);
-	calculate_derivs(i,   j+1, k+1, ivar, nvar, n1, n2, n3, v, vv, 7);
-	calculate_derivs(i+1, j+1, k+1, ivar, nvar, n1, n2, n3, v, vv, 8);
-	v_intpol[2]=interpol(ha*(a-1),hb*b,    hp*c,    vv,2);
-	v_intpol[3]=interpol(ha*a,    hb*(b-1),hp*c,    vv,3);
-	v_intpol[4]=interpol(ha*(a-1),hb*(b-1),hp*c,    vv,4);
-	v_intpol[5]=interpol(ha*a,    hb*b,    hp*(c-1),vv,5);
-	v_intpol[6]=interpol(ha*(a-1),hb*b,    hp*(c-1),vv,6);
-	v_intpol[7]=interpol(ha*a,    hb*(b-1),hp*(c-1),vv,7);
-	v_intpol[8]=interpol(ha*(a-1),hb*(b-1),hp*(c-1),vv,8);
-	for(j=1;j<=8;j++)
-		printf("Value U[V] at origin:%16.15f\tj=%d\n",-2*v_intpol[j],j);*/
-/*	result = 0;
-	for(mi=i-npol; mi<=i+npol; mi++){
-		al_m  = Pih*(2*mi+1)/n1;
-		for(mj=j-npol; mj<=j+npol; mj++){
-			be_m  = Pih*(2*mj+1)/n2;  
-			for(mk=k-npol; mk<=k+npol; mk++){
-				phi_m = 2.*Pi*mk/n3;
-				g_m   = 1;
-				for(ni=i-npol; ni<=i+npol; ni++){
-					al_n  = Pih*(2*ni+1)/n1;
-					if(ni != mi) g_m *= (al - al_n)/(al_m - al_n);
-				}
-				for(nj=j-npol; nj<=j+npol; nj++){
-					be_n  = Pih*(2*nj+1)/n2;  
-					if(nj != mj) g_m *= (be - be_n)/(be_m - be_n);
-				}
-				for(nk=k-npol; nk<=k+npol; nk++){
-					phi_n = 2.*Pi*nk/n3;
-					if(nk != mk) g_m *= (phi-phi_n)/(phi_m - phi_n);
-				}
-				result += g_m*v.d0[Index(ivar,mi,mj,mk,nvar,n1,n2,n3)];
-			}
-		}
-	}	*/
+// -------------------------------------------------------------------------------
+// Calculates the value of v at an arbitrary position (x,y,z)
+double
+PunctIntPolAtArbitPosition (int ivar, int nvar, int n1,
+			    int n2, int n3, derivs v, double x, double y,
+			    double z)
+{
+  DECLARE_CCTK_PARAMETERS;
+  double xs, ys, zs, rs2, phi, X, R, A, B, aux1, aux2, result, Ui;
+
+  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 (sqrt (-aux1 + aux2));
+  if (x < 0)
+    R = Pi - R;
+
+  A = 2 * tanh (0.5 * X) - 1;
+  B = tan (0.5 * R - Piq);
+
+  result = PunctEvalAtArbitPosition (v.d0, A, B, phi, n1, n2, n3);
+
+  Ui = (A - 1) * result;
+
+  return Ui;
 }
 
 // -------------------------------------------------------------------------------
diff --git a/src/TwoPunctures.c b/src/TwoPunctures.c
index b835145..519aabb 100644
--- a/src/TwoPunctures.c
+++ b/src/TwoPunctures.c
@@ -32,17 +32,20 @@ TwoPunctures (CCTK_ARGUMENTS)
   int nvar = 1, n1 = npoints_A, n2 = npoints_B, n3 = npoints_phi;
 
   int i, j, k, ntotal = n1 * n2 * n3 * nvar;
-  double *F;
-  derivs u, v;
+  static double *F = NULL;
+  static derivs u, v;
 
-  F = dvector (0, ntotal - 1);
-  allocate_derivs (&u, ntotal);
-  allocate_derivs (&v, ntotal);
+  if (! F) {
+    /* Solve only when called for the first time */
+    F = dvector (0, ntotal - 1);
+    allocate_derivs (&u, ntotal);
+    allocate_derivs (&v, ntotal);
 
-  CCTK_INFO ("Solving puncture equation");
-  Newton (nvar, n1, n2, n3, v, 1.e-10, 5);
+    CCTK_INFO ("Solving puncture equation");
+    Newton (nvar, n1, n2, n3, v, Newton_tol, Newton_maxit);
 
-  F_of_v (nvar, n1, n2, n3, v, F, u);
+    F_of_v (nvar, n1, n2, n3, v, F, u);
+  }
 
   CCTK_INFO ("Interpolating result");
   if (CCTK_EQUALS(metric_type, "static conformal")) {
@@ -71,8 +74,10 @@ TwoPunctures (CCTK_ARGUMENTS)
         const double r_minus
           = sqrt(pow2(x[ind] + par_b) + pow2(y[ind]) + pow2(z[ind]));
         
-        const double U = PunctIntPolAtArbitPosition
+        const double U = PunctTaylorExpandAtArbitPosition
 	  (0, nvar, n1, n2, n3, v, x[ind], y[ind], z[ind]);
+/*         const double U = PunctIntPolAtArbitPosition */
+/* 	  (0, nvar, n1, n2, n3, v, x[ind], y[ind], z[ind]); */
         const double psi1 = 1
           + 0.5 * par_m_plus / r_plus
           + 0.5 * par_m_minus / r_minus + U;
@@ -158,7 +163,7 @@ TwoPunctures (CCTK_ARGUMENTS)
             psizz[ind] = pzz / static_psi;
           }
 
-        } /* if conformal-state>0 */
+        } /* if conformal-state > 0 */
           
         puncture_u[ind] = U;
 
@@ -180,7 +185,10 @@ TwoPunctures (CCTK_ARGUMENTS)
     }
   }
 
-  free_dvector (F, 0, ntotal - 1);
-  free_derivs (&u, ntotal);
-  free_derivs (&v, ntotal);
+  if (0) {
+    /* Keep the result around for the next time */
+    free_dvector (F, 0, ntotal - 1);
+    free_derivs (&u, ntotal);
+    free_derivs (&v, ntotal);
+  }
 }
diff --git a/src/TwoPunctures.h b/src/TwoPunctures.h
index ab0d841..9915181 100644
--- a/src/TwoPunctures.h
+++ b/src/TwoPunctures.h
@@ -43,6 +43,9 @@ double PunctEvalAtArbitPosition (double *v, double A, double B, double phi,
 void calculate_derivs (int i, int j, int k, int ivar, int nvar, int n1,
 		       int n2, int n3, derivs v, derivs vv);
 double interpol (double a, double b, double c, derivs v);
+double PunctTaylorExpandAtArbitPosition (int ivar, int nvar, int n1,
+                                         int n2, int n3, derivs v, double x,
+                                         double y, double z);
 double PunctIntPolAtArbitPosition (int ivar, int nvar, int n1,
 				   int n2, int n3, derivs v, double x,
 				   double y, double z);
diff --git a/test/bam.par b/test/bam.par
new file mode 100644
index 0000000..6c401c4
--- /dev/null
+++ b/test/bam.par
@@ -0,0 +1,61 @@
+# $Header$
+
+ActiveThorns = "Boundary CartGrid3D Time CoordBase SymBase"
+ActiveThorns = "PUGH PUGHReduce PUGHSlab"
+ActiveThorns = "IOASCII IOBasic IOUtil"
+ActiveThorns = "ADMBase StaticConformal"
+ActiveThorns = "ADMConstraints ADMCoupling ADMMacros SpaceMask"
+ActiveThorns = "BAM_Elliptic EllBase"
+
+Cactus::cctk_itlast		= 0
+ 
+Time::dtfac			= 0.25
+
+Grid::type			= "byrange"
+Grid::domain			= "full"
+Grid::xmin                      = -4
+Grid::xmax                      =  4
+Grid::ymin                      = -4
+Grid::ymax			=  4
+Grid::zmin                      = -4
+Grid::zmax			=  4
+Driver::global_nx		= 35
+Driver::global_ny		= 35
+Driver::global_nz		= 35
+
+ADMBase::initial_data		= "BrBr"
+
+BAM_Elliptic::redblack_decouple = "yes"
+BAM_Elliptic::bam_persist       = "yes"
+BAM_Elliptic::brbr_tolres2  	= 1e-6
+BAM_Elliptic::brbr_maxcycles 	= 100
+BAM_Elliptic::brbr_mincycles 	= 3
+
+BAM_Elliptic::bhm1		=  1.5
+BAM_Elliptic::bhx1		=  1.5
+BAM_Elliptic::bhpy1		=  2.0
+BAM_Elliptic::bhsy1		=  0.5
+BAM_Elliptic::bhsz1		= -0.5
+
+BAM_Elliptic::bhm2		=  1.0
+BAM_Elliptic::bhx2		= -1.5
+BAM_Elliptic::bhpy2		= -2.0
+BAM_Elliptic::bhsx2		= -1.0
+BAM_Elliptic::bhsz2		= -1.0
+
+BAM_Elliptic::bam_bound		= "bamrobin"
+
+ADMBase::metric_type		= "static conformal"
+StaticConformal::conformal_storage = "factor+derivs+2nd derivs"
+
+ADMBase::lapse_evolution_method	= "static"
+ADMBase::initial_lapse		= "one"
+  
+IO::out_dir			= $parfile
+IO::out_fileinfo		= "axis labels"
+
+IOBasic::outScalar_every	= 1
+IOBasic::outScalar_vars		= "staticconformal::psi admbase::metric admbase::curv admconstraints::hamiltonian admconstraints::momentum"
+
+IOASCII::out1D_every		= 1
+IOASCII::out1D_vars		= "staticconformal::psi admbase::metric admbase::curv admconstraints::hamiltonian admconstraints::momentum"
diff --git a/test/idpuncture.par b/test/idpuncture.par
new file mode 100644
index 0000000..96df1aa
--- /dev/null
+++ b/test/idpuncture.par
@@ -0,0 +1,64 @@
+# $Header$
+
+ActiveThorns = "Boundary CartGrid3D Time CoordBase SymBase"
+ActiveThorns = "PUGH PUGHReduce PUGHSlab"
+ActiveThorns = "IOASCII IOBasic IOUtil"
+ActiveThorns = "ADMBase StaticConformal"
+ActiveThorns = "ADMConstraints ADMCoupling ADMMacros SpaceMask"
+ActiveThorns = "IDPuncture TATelliptic TATJacobi TATMG TATPETSc"
+
+Cactus::cctk_itlast		= 0
+ 
+Time::dtfac			= 0.25
+
+Grid::type			= "byrange"
+Grid::domain			= "full"
+Grid::xmin                      = -4
+Grid::xmax                      =  4
+Grid::ymin                      = -4
+Grid::ymax			=  4
+Grid::zmin                      = -4
+Grid::zmax			=  4
+Driver::global_nx		= 35
+Driver::global_ny		= 35
+Driver::global_nz		= 35
+
+ADMBase::initial_data		= "puncture"
+
+IDPuncture::verbose		= "yes"
+
+IDPuncture::bhm1		=  1.5
+IDPuncture::bhx1		=  1.5
+IDPuncture::bhpy1		=  2.0
+IDPuncture::bhsy1		=  0.5
+IDPuncture::bhsz1		= -0.5
+
+IDPuncture::bhm2		=  1.0
+IDPuncture::bhx2		= -1.5
+IDPuncture::bhpy2		= -2.0
+IDPuncture::bhsx2		= -1.0
+IDPuncture::bhsz2		= -1.0
+
+IDPuncture::bound		= "robin"
+
+IDPuncture::solver		= "TATPETSc"
+
+ADMBase::metric_type		= "static conformal"
+StaticConformal::conformal_storage = "factor+derivs+2nd derivs"
+
+ADMBase::lapse_evolution_method	= "static"
+ADMBase::initial_lapse		= "one"
+
+TATJacobi::verbose		= yes
+
+TATPETSc::verbose		= yes
+TATPETsc::options		= "-snes_monitor -get_total_flops"
+  
+IO::out_dir			= $parfile
+IO::out_fileinfo		= "axis labels"
+
+IOBasic::outScalar_every	= 1
+IOBasic::outScalar_vars		= "staticconformal::psi admbase::metric admbase::curv admconstraints::hamiltonian admconstraints::momentum idpuncture::u idpuncture::M idpuncture::N"
+
+IOASCII::out1D_every		= 1
+IOASCII::out1D_vars		= "staticconformal::psi admbase::metric admbase::curv admconstraints::hamiltonian admconstraints::momentum idpuncture::u idpuncture::M idpuncture::N"
diff --git a/test/twopunctures.par b/test/twopunctures.par
index 85be41d..29af25b 100644
--- a/test/twopunctures.par
+++ b/test/twopunctures.par
@@ -1,6 +1,11 @@
 # $Header$
 
-ActiveThorns = "Boundary CartGrid3D Time CoordBase SymBase PUGH PUGHReduce PUGHSlab IOASCII IOBasic IOUtil ADMBase StaticConformal TwoPunctures"
+ActiveThorns = "Boundary CartGrid3D Time CoordBase SymBase"
+ActiveThorns = "PUGH PUGHReduce PUGHSlab"
+ActiveThorns = "IOASCII IOBasic IOUtil"
+ActiveThorns = "ADMBase StaticConformal"
+ActiveThorns = "ADMConstraints ADMCoupling ADMMacros SpaceMask"
+ActiveThorns = "TwoPunctures"
 
 Cactus::cctk_itlast		= 0
 
@@ -26,17 +31,17 @@ TwoPunctures::npoints_A         = 20
 TwoPunctures::npoints_B         = 20
 TwoPunctures::npoints_phi       = 12
 
-TwoPunctures::par_b             = 1.5
+TwoPunctures::par_b             =  1.5
 
-TwoPunctures::par_m_plus        = 1.5
-TwoPunctures::par_P_plus[1]     = 2.0
-TwoPunctures::par_S_plus[1]     = 0.5
-TwoPunctures::par_S_plus[2]     = 0.5
+TwoPunctures::par_m_plus        =  1.5
+TwoPunctures::par_P_plus[1]     =  2.0
+TwoPunctures::par_S_plus[1]     =  0.5
+TwoPunctures::par_S_plus[2]     = -0.5
 
 TwoPunctures::par_m_minus       =  1.0
 TwoPunctures::par_P_minus[1]    = -2.0
 TwoPunctures::par_S_minus[0]    = -1.0
-TwoPunctures::par_S_minus[2]    =  1.0
+TwoPunctures::par_S_minus[2]    = -1.0
 
 ADMBase::metric_type		= "static conformal"
 StaticConformal::conformal_storage = "factor+derivs+2nd derivs"
@@ -48,7 +53,7 @@ IO::out_dir			= $parfile
 IO::out_fileinfo		= "axis labels"
 
 IOBasic::outScalar_every	= 1
-IOBasic::outScalar_vars		= "staticconformal::psi admbase::metric admbase::curv TwoPunctures::puncture_u"
+IOBasic::outScalar_vars		= "staticconformal::psi admbase::metric admbase::curv admconstraints::hamiltonian admconstraints::momentum TwoPunctures::puncture_u"
 
 IOASCII::out1D_every		= 1
-IOASCII::out1D_vars		= "staticconformal::psi admbase::metric admbase::curv TwoPunctures::puncture_u"
+IOASCII::out1D_vars		= "staticconformal::psi admbase::metric admbase::curv admconstraints::hamiltonian admconstraints::momentum TwoPunctures::puncture_u"