From c5078deec69fa8bba4d130d3fadb2d0d083b21e2 Mon Sep 17 00:00:00 2001 From: jese Date: Wed, 26 Mar 2008 15:48:41 +0000 Subject: Version which works for runs and calculates bare_mass, but may not properly set calculated bare masses. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/TwoPunctures/trunk@83 b2a53a04-0f4f-0410-87ed-f9f25ced00cf --- src/TwoPunctures.c | 66 +++++++++++++++++++++++++----------------------------- 1 file changed, 31 insertions(+), 35 deletions(-) diff --git a/src/TwoPunctures.c b/src/TwoPunctures.c index 8ef7911..bc70973 100644 --- a/src/TwoPunctures.c +++ b/src/TwoPunctures.c @@ -243,19 +243,22 @@ TwoPunctures (CCTK_ARGUMENTS) { set_initial_guess(cctkGH, v); } + mp = *(CCTK_REAL *) CCTK_ParameterGet("par_m_plus", "twopunctures", NULL); + mm = *(CCTK_REAL *) CCTK_ParameterGet("par_m_minus", "twopunctures", NULL); - if(!give_bare_mass) { - Newton (cctkGH, nvar, n1, n2, n3, v, Newton_tol, 2); + if(!(give_bare_mass)) { + Newton (cctkGH, nvar, n1, n2, n3, v, Newton_tol, 2); - F_of_v (cctkGH, nvar, n1, n2, n3, v, F, u); + F_of_v (cctkGH, nvar, n1, n2, n3, v, F, u); - - up = PunctTaylorExpandAtArbitPosition + up = PunctIntPolAtArbitPosition (0, nvar, n1, n2, n3, v, par_b, 0., 0.); - um = PunctTaylorExpandAtArbitPosition + um = PunctIntPolAtArbitPosition (0, nvar, n1, n2, n3, v, -par_b, 0., 0.); - mp = target_M_plus; - mm = target_M_minus; + + mp = *(CCTK_REAL *) CCTK_ParameterGet("target_M_plus", "twopunctures", NULL); + mm = *(CCTK_REAL *) CCTK_ParameterGet("target_M_minus", "twopunctures", NULL); + M_p= target_M_plus; M_m= target_M_minus; for(l=0; l<32;l++) { @@ -268,15 +271,11 @@ TwoPunctures (CCTK_ARGUMENTS) } char valbuf_p[100], valbuf_m[100]; - snprintf (valbuf_p,sizeof valbuf_p, "%g", (float) mp); + sprintf (valbuf_p,"%f", (float) mp); CCTK_ParameterSet ("par_m_plus", "twopunctures", valbuf_p); - snprintf (valbuf_m,sizeof valbuf_m, "%g", (float) mm); + sprintf (valbuf_m,"%f", (float) mm); CCTK_ParameterSet ("par_m_minus", "twopunctures", valbuf_m); } - else { - mp = par_m_plus; - mm = par_m_minus; - } Newton (cctkGH, nvar, n1, n2, n3, v, Newton_tol, Newton_maxit); @@ -286,18 +285,15 @@ TwoPunctures (CCTK_ARGUMENTS) "The two puncture masses are %g and %g", (double) mm, (double) mp); - + CCTK_VInfo (CCTK_THORNSTRING, + "The puncture masses are %g and %g", + (double) par_m_minus, (double) par_m_plus); + /* 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, 1, 0, 0, n1, n2, n3)); CCTK_VInfo (CCTK_THORNSTRING, "ADM mass is %g\n", (double)admMass); } - else { - - mp = par_m_plus; - mm = par_m_minus; - - } if (CCTK_EQUALS(grid_setup_method, "Taylor expansion")) { @@ -416,21 +412,21 @@ TwoPunctures (CCTK_ARGUMENTS) if (r_minus < TP_Tiny) r_minus = TP_Tiny; CCTK_REAL psi1 = 1 - + 0.5 * mp / r_plus - + 0.5 * mm / r_minus + U; + + 0.5 * par_m_plus / r_plus + + 0.5 * par_m_minus / r_minus + U; #define EXTEND(M,r) \ ( M * (3./8 * pow(r, 4) / pow(TP_Extend_Radius, 5) - \ 5./4 * pow(r, 2) / pow(TP_Extend_Radius, 3) + \ 15./8 / TP_Extend_Radius)) if (r_plus < TP_Extend_Radius) { psi1 = 1 - + 0.5 * EXTEND(mp,r_plus) - + 0.5 * mm / r_minus + U; + + 0.5 * EXTEND(par_m_plus,r_plus) + + 0.5 * par_m_minus / r_minus + U; } if (r_minus < TP_Extend_Radius) { psi1 = 1 - + 0.5 * EXTEND(mm,r_minus) - + 0.5 * mp / r_plus + U; + + 0.5 * EXTEND(par_m_minus,r_minus) + + 0.5 * par_m_plus / r_plus + U; } CCTK_REAL static_psi = 1; @@ -464,7 +460,7 @@ TwoPunctures (CCTK_ARGUMENTS) ir = EXTEND(1., rp); } - s1 = 0.5*mp*ir; + s1 = 0.5*par_m_plus*ir; s3 = -s1*ir*ir; s5 = -3.0*s3*ir*ir; @@ -495,7 +491,7 @@ TwoPunctures (CCTK_ARGUMENTS) ir = EXTEND(1., rp); } - s1 = 0.5*mm*ir; + s1 = 0.5*par_m_minus*ir; s3 = -s1*ir*ir; s5 = -3.0*s3*ir*ir; @@ -555,18 +551,18 @@ TwoPunctures (CCTK_ARGUMENTS) if (antisymmetric_lapse || averaged_lapse) { alp[ind] = - ((1.0 -0.5*mp/r_plus -0.5*mm/r_minus) - /(1.0 +0.5*mp/r_plus +0.5*mm/r_minus)); + ((1.0 -0.5*par_m_plus/r_plus -0.5*par_m_minus/r_minus) + /(1.0 +0.5*par_m_plus/r_plus +0.5*par_m_minus/r_minus)); if (r_plus < TP_Extend_Radius) { alp[ind] = - ((1.0 -0.5*EXTEND(mp, r_plus) -0.5*mm/r_minus) - /(1.0 +0.5*EXTEND(mp, r_plus) +0.5*mm/r_minus)); + ((1.0 -0.5*EXTEND(par_m_plus, r_plus) -0.5*par_m_minus/r_minus) + /(1.0 +0.5*EXTEND(par_m_plus, r_plus) +0.5*par_m_minus/r_minus)); } if (r_minus < TP_Extend_Radius) { alp[ind] = - ((1.0 -0.5*EXTEND(mm, r_minus) -0.5*mp/r_plus) - /(1.0 +0.5*EXTEND(mm, r_minus) +0.5*mp/r_plus)); + ((1.0 -0.5*EXTEND(par_m_minus, r_minus) -0.5*par_m_plus/r_plus) + /(1.0 +0.5*EXTEND(par_m_minus, r_minus) +0.5*par_m_plus/r_plus)); } if (averaged_lapse) { -- cgit v1.2.3