[Commits] [svn:einsteintoolkit] TwoPunctures/trunk/ (Rev. 124)

User: hinder Date: 2011/09/25 04:43 PM

Added: /trunk/src/ Metadata.cc

Modified: /trunk/ interface.ccl, param.ccl, schedule.ccl /trunk/src/ TwoPunctures.c, make.code.defn

Log: 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.

File Changes:

Directory: /trunk/src/ ======================

File [added]: Metadata.cc Delta lines: +64 -0 =================================================================== --- trunk/src/Metadata.cc (rev 0) +++ trunk/src/Metadata.cc 2011-09-25 21:43:44 UTC (rev 124) @@ -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; + } +}

File [modified]: TwoPunctures.c Delta lines: +37 -13 =================================================================== --- trunk/src/TwoPunctures.c 2011-05-20 11:02:32 UTC (rev 123) +++ trunk/src/TwoPunctures.c 2011-09-25 21:43:44 UTC (rev 124) @@ -220,7 +220,7 @@ 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 @@ 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 @@ 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 @@ 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 @@ 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"))

File [modified]: make.code.defn Delta lines: +1 -1 =================================================================== --- trunk/src/make.code.defn 2011-05-20 11:02:32 UTC (rev 123) +++ trunk/src/make.code.defn 2011-09-25 21:43:44 UTC (rev 124) @@ -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 =

Directory: /trunk/ ==================

File [modified]: interface.ccl Delta lines: +16 -0 =================================================================== --- trunk/interface.ccl 2011-05-20 11:02:32 UTC (rev 123) +++ trunk/interface.ccl 2011-09-25 21:43:44 UTC (rev 124) @@ -40,6 +40,18 @@ 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 @@ 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)"

File [modified]: param.ccl Delta lines: +3 -0 =================================================================== --- trunk/param.ccl 2011-05-20 11:02:32 UTC (rev 123) +++ trunk/param.ccl 2011-09-25 21:43:44 UTC (rev 124) @@ -23,8 +23,11 @@

USES KEYWORD conformal_storage

+SHARES: IO

+USES STRING out_dir

+ RESTRICTED:

BOOLEAN verbose "Print screen output while solving"

File [modified]: schedule.ccl Delta lines: +21 -12 =================================================================== --- trunk/schedule.ccl 2011-05-20 11:02:32 UTC (rev 123) +++ trunk/schedule.ccl 2011-09-25 21:43:44 UTC (rev 124) @@ -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 (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" }