From 32ae652a29c362e93582fbbb718a37152a09cf33 Mon Sep 17 00:00:00 2001 From: hinder Date: Sun, 25 Sep 2011 21:43:44 +0000 Subject: 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 --- interface.ccl | 16 ++++++++++++++ param.ccl | 3 +++ schedule.ccl | 33 ++++++++++++++++++---------- src/Metadata.cc | 64 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ src/TwoPunctures.c | 50 +++++++++++++++++++++++++++++++----------- src/make.code.defn | 2 +- 6 files changed, 142 insertions(+), 26 deletions(-) create mode 100644 src/Metadata.cc 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 +#include +#include + +#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 + + 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 = -- cgit v1.2.3