Add output of metadata file TwoPunctures.bbh

The format and keys are defined in http://arxiv.org/abs/0709.0093 with additional draft keys defined for the NR-AR project (to be committed to the arXiv version when finalized).

Also introduced grid scalars for the ADM energy and angular momentum and the puncture ADM masses.

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

6 files changed, 142 insertions, 26 deletions

diff --git a/interface.ccl b/interface.ccl
index e58faac..668cfc0 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -40,6 +40,18 @@ USES FUNCTION Set_Rho_ADM
 USES FUNCTION Set_Initial_Guess_for_u
 USES FUNCTION Rescale_Sources
 
+PRIVATE:
+
+CCTK_REAL energy TYPE=scalar
+{
+  E
+} "ADM energy of the Bowen-York spacetime"
+
+CCTK_REAL angular_momentum TYPE=scalar
+{
+  J1, J2, J3
+} "Angular momentum of the Bowen-York spacetime"
+
 PUBLIC:
 
 CCTK_REAL bare_mass TYPE=scalar
@@ -47,4 +59,8 @@ CCTK_REAL bare_mass TYPE=scalar
   mp mm
 } "Bare masses of the punctures"
 
+CCTK_REAL puncture_adm_mass TYPE=scalar
+{
+  mp_adm mm_adm
+} "ADM masses of the punctures (measured at the other spatial infinities)"
 
diff --git a/param.ccl b/param.ccl
index 1978c35..b4108aa 100644
--- a/param.ccl
+++ b/param.ccl
@@ -23,6 +23,9 @@ SHARES: StaticConformal
 
 USES KEYWORD conformal_storage
 
+SHARES: IO
+
+USES STRING out_dir
 
 
 RESTRICTED:
diff --git a/schedule.ccl b/schedule.ccl
index 4353c66..643ddb5 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -2,6 +2,8 @@
 
 if (CCTK_Equals(initial_data, "twopunctures"))
 {
+  STORAGE: energy, angular_momentum, puncture_adm_mass
+
   if (keep_u_around) {
     STORAGE: puncture_u bare_mass
   }
@@ -13,20 +15,27 @@ if (CCTK_Equals(initial_data, "twopunctures"))
 
   if (schedule_in_ADMBase_InitialData)
   {
-    SCHEDULE TwoPunctures IN ADMBase_InitialData
-    {
-      LANG: C
-      STORAGE: puncture_u bare_mass
-      # SYNC: ADMBase::metric ADMBase::curv ADMBase::lapse
-    } "Create puncture black hole initial data"
+      SCHEDULE GROUP TwoPunctures_Group IN ADMBase_InitialData
+      {
+      } "TwoPunctures initial data group"
   }
   else
   {
-    SCHEDULE TwoPunctures AT Initial AFTER ADMBase_InitialData BEFORE ADMBase_PostInitial AFTER HydroBase_Initial before SetTmunu before HydroBase_Prim2ConInitial
-    {
-      LANG: C
-      STORAGE: puncture_u bare_mass
-      # SYNC: ADMBase::metric ADMBase::curv ADMBase::lapse
-    } "Create puncture black hole initial data"
+      SCHEDULE GROUP TwoPunctures_Group AT Initial AFTER ADMBase_InitialData BEFORE ADMBase_PostInitial AFTER HydroBase_Initial before SetTmunu before HydroBase_Prim2ConInitial
+      {
+      } "TwoPunctures initial data group"
   }
+
+  SCHEDULE TwoPunctures IN TwoPunctures_Group
+  {
+    LANG: C
+    STORAGE: puncture_u bare_mass
+    # SYNC: ADMBase::metric ADMBase::curv ADMBase::lapse
+  } "Create puncture black hole initial data"
+
+  SCHEDULE TwoPunctures_Metadata IN TwoPunctures_Group after TwoPunctures
+  {
+    LANG: C
+    OPTIONS: global
+  } "Output TwoPunctures metadata"
 }
diff --git a/src/Metadata.cc b/src/Metadata.cc
new file mode 100644
index 0000000..b67903b
--- /dev/null
+++ b/src/Metadata.cc
@@ -0,0 +1,64 @@
+
+#include <iostream>
+#include <fstream>
+#include <iomanip>
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+using namespace std;
+
+extern "C"
+void TwoPunctures_Metadata (CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+
+  if (CCTK_MyProc(cctkGH) == 0)
+  {
+    ofstream o;
+    o.open(string(string(out_dir) + "/TwoPunctures.bbh").c_str());
+
+    o << setprecision(19);
+
+    o << "\
+# ==================================\n\
+# Numerical Relativity Metadata file\n\
+# ==================================\n\
+#\n\
+# This file contains information about the simulation provided by the\n\
+# TwoPunctures thorn.  The format is described in the NR Data Format Document\n\
+# http://arxiv.org/abs/0709.0093 [draft SVN r707].\n\
+" << endl;
+
+    o << "[metadata]" << endl;
+    o << "initial-ADM-energy            = " << *E << endl;
+    o << "initial-ADM-angular-momentumx = " << *J1 << endl;
+    o << "initial-ADM-angular-momentumy = " << *J2 << endl;
+    o << "initial-ADM-angular-momentumz = " << *J3 << endl;
+    o << "initial-separation            = " << par_b * 2 << endl;
+    o << "initial-data-type             = Bowen-York" << endl;
+    o << "initial-data-bibtex-keys      = Bowen:1980yu Brandt:1997tf Ansorg:2004ds" << endl;
+    o << "initial-bh-position1x         = " << par_b + center_offset[0] << endl;
+    o << "initial-bh-position1y         = " << center_offset[1] << endl;
+    o << "initial-bh-position1z         = " << center_offset[2] << endl;
+    o << "initial-bh-position2x         = " << -par_b + center_offset[0] << endl;
+    o << "initial-bh-position2y         = " << center_offset[1] << endl;
+    o << "initial-bh-position2z         = " << center_offset[2] << endl;
+    o << "initial-bh-momentum1x         = " << par_P_plus[0] << endl;
+    o << "initial-bh-momentum1y         = " << par_P_plus[1] << endl;
+    o << "initial-bh-momentum1z         = " << par_P_plus[2] << endl;
+    o << "initial-bh-momentum2x         = " << par_P_minus[0] << endl;
+    o << "initial-bh-momentum2y         = " << par_P_minus[1] << endl;
+    o << "initial-bh-momentum2z         = " << par_P_minus[2] << endl;
+    o << "initial-bh-spin1x             = " << par_S_plus[0] << endl;
+    o << "initial-bh-spin1y             = " << par_S_plus[1] << endl;
+    o << "initial-bh-spin1z             = " << par_S_plus[2] << endl;
+    o << "initial-bh-spin2x             = " << par_S_minus[0] << endl;
+    o << "initial-bh-spin2y             = " << par_S_minus[1] << endl;
+    o << "initial-bh-spin2z             = " << par_S_minus[2] << endl;
+    o << "initial-bh-puncture-adm-mass1 = " << *mp_adm << endl;
+    o << "initial-bh-puncture-adm-mass2 = " << *mm_adm << endl;
+  }
+}
diff --git a/src/TwoPunctures.c b/src/TwoPunctures.c
index f55380f..8835c74 100644
--- a/src/TwoPunctures.c
+++ b/src/TwoPunctures.c
@@ -220,7 +220,7 @@ TwoPunctures (CCTK_ARGUMENTS)
   CCTK_REAL admMass;
 
   if (! F) {
-    CCTK_REAL Mp_adm, Mm_adm, up, um;
+    CCTK_REAL up, um;
     /* Solve only when called for the first time */
     F = dvector (0, ntotal - 1);
     allocate_derivs (&u, ntotal);
@@ -257,7 +257,7 @@ TwoPunctures (CCTK_ARGUMENTS)
        target ADM masses target_M_plus and target_M_minus and with initial 
        guesses given by par_m_plus and par_m_minus. */
     if(!(give_bare_mass)) {
-      CCTK_REAL tmp, Mp_adm_err, Mm_adm_err;
+      CCTK_REAL tmp, mp_adm_err, mm_adm_err;
       char valbuf[100];
 
       CCTK_REAL M_p = target_M_plus;
@@ -280,14 +280,14 @@ TwoPunctures (CCTK_ARGUMENTS)
         um = PunctIntPolAtArbitPosition(0, nvar, n1, n2, n3, v,-par_b, 0., 0.);
 
         /* Calculate the ADM masses from the current bare mass guess */
-        Mp_adm = (1 + up) * *mp + *mp * *mm / (4. * par_b);
-        Mm_adm = (1 + um) * *mm + *mp * *mm / (4. * par_b);
+        *mp_adm = (1 + up) * *mp + *mp * *mm / (4. * par_b);
+        *mm_adm = (1 + um) * *mm + *mp * *mm / (4. * par_b);
 
         /* Check how far the current ADM masses are from the target */
-        Mp_adm_err = fabs(M_p-Mp_adm);
-        Mm_adm_err = fabs(M_m-Mm_adm);
+        mp_adm_err = fabs(M_p-*mp_adm);
+        mm_adm_err = fabs(M_m-*mm_adm);
         CCTK_VInfo (CCTK_THORNSTRING, "ADM mass error: M_p_err=%.4g, M_m_err=%.4g",
-                    (double)Mp_adm_err, (double)Mm_adm_err);
+                    (double) mp_adm_err, (double) mm_adm_err);
         
         /* Invert the ADM mass equation and update the bare mass guess so that
            it gives the correct target ADM masses */
@@ -304,8 +304,8 @@ TwoPunctures (CCTK_ARGUMENTS)
         sprintf (valbuf,"%.17g", (double) *mm);
         CCTK_ParameterSet ("par_m_minus", "twopunctures", valbuf);
         
-      } while ( (Mp_adm_err > adm_tol) ||
-                (Mm_adm_err > adm_tol) );
+      } while ( (mp_adm_err > adm_tol) ||
+                (mm_adm_err > adm_tol) );
                 
       CCTK_VInfo (CCTK_THORNSTRING, "Found bare masses.");
     }
@@ -322,16 +322,40 @@ TwoPunctures (CCTK_ARGUMENTS)
     um = PunctIntPolAtArbitPosition(0, nvar, n1, n2, n3, v,-par_b, 0., 0.);
 
     /* Calculate the ADM masses from the current bare mass guess */
-    Mp_adm = (1 + up) * *mp + *mp * *mm / (4. * par_b);
-    Mm_adm = (1 + um) * *mm + *mp * *mm / (4. * par_b);
+    *mp_adm = (1 + up) * *mp + *mp * *mm / (4. * par_b);
+    *mm_adm = (1 + um) * *mm + *mp * *mm / (4. * par_b);
 
-    CCTK_VInfo (CCTK_THORNSTRING, "Puncture 1 ADM mass is %g", (double)Mp_adm);
-    CCTK_VInfo (CCTK_THORNSTRING, "Puncture 2 ADM mass is %g", (double)Mm_adm);
+    CCTK_VInfo (CCTK_THORNSTRING, "Puncture 1 ADM mass is %g", (double) *mp_adm);
+    CCTK_VInfo (CCTK_THORNSTRING, "Puncture 2 ADM mass is %g", (double) *mm_adm);
 
     /* print out ADM mass, eq.: \Delta M_ADM=2*r*u=4*b*V for A=1,B=0,phi=0 */
     admMass = (*mp + *mm
                - 4*par_b*PunctEvalAtArbitPosition(v.d0, 0, 1, 0, 0, nvar, n1, n2, n3));
     CCTK_VInfo (CCTK_THORNSTRING, "The total ADM mass is %g", (double) admMass);
+    *E = admMass;
+
+    /*
+      Run this in Mathematica (version 8 or later) with
+        math -script <file>
+
+      Needs["SymbolicC`"];
+      co = Table["center_offset[" <> ToString[i] <> "]", {i, 0, 2}];
+      r1 = co + {"par_b", 0, 0};
+      r2 = co + {-"par_b", 0, 0};
+      {p1, p2} = Table["par_P_" <> bh <> "[" <> ToString[i] <> "]", {bh, {"plus", "minus"}}, {i, 0, 2}];
+      {s1, s2} = Table["par_S_" <> bh <> "[" <> ToString[i] <> "]", {bh, {"plus", "minus"}}, {i, 0, 2}];
+
+      J = Cross[r1, p1] + Cross[r2, p2] + s1 + s2;
+
+      JVar = Table["*J" <> ToString[i], {i, 1, 3}];
+      Print[OutputForm@StringReplace[
+        ToCCodeString@MapThread[CAssign[#1, CExpression[#2]] &, {JVar, J}], 
+        "\"" -> ""]];
+     */
+
+    *J1 = -(center_offset[2]*par_P_minus[1]) + center_offset[1]*par_P_minus[2] - center_offset[2]*par_P_plus[1] + center_offset[1]*par_P_plus[2] + par_S_minus[0] + par_S_plus[0];
+    *J2 = center_offset[2]*par_P_minus[0] - center_offset[0]*par_P_minus[2] + par_b*par_P_minus[2] + center_offset[2]*par_P_plus[0] - center_offset[0]*par_P_plus[2] - par_b*par_P_plus[2] + par_S_minus[1] + par_S_plus[1];
+    *J3 = -(center_offset[1]*par_P_minus[0]) + center_offset[0]*par_P_minus[1] - par_b*par_P_minus[1] - center_offset[1]*par_P_plus[0] + center_offset[0]*par_P_plus[1] + par_b*par_P_plus[1] + par_S_minus[2] + par_S_plus[2];
   }
 
   if (CCTK_EQUALS(grid_setup_method, "Taylor expansion"))
diff --git a/src/make.code.defn b/src/make.code.defn
index b0ebf81..14e70ea 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -1,7 +1,7 @@
 # Main make.code.defn file for thorn TwoPunctures
 
 # Source files in this directory
-SRCS = CoordTransf.c Equations.c FuncAndJacobian.c Newton.c TwoPunctures.c TP_utilities.c ParamCheck.c
+SRCS = CoordTransf.c Equations.c FuncAndJacobian.c Newton.c TwoPunctures.c TP_utilities.c ParamCheck.c Metadata.cc
 
 # Subdirectories containing source files
 SUBDIRS =