[Commits] [svn:einsteintoolkit] incoming/TOVSolverHot/ (Rev. 25)

User: cott Date: 2010/12/04 03:53 PM

Modified: /TOVSolverHot/ param.ccl /TOVSolverHot/src/ external.inc, tov.c

Log: * there may be only one K and Gamma; Make parameters scalar.

File Changes:

Directory: /TOVSolverHot/ =========================

File [modified]: param.ccl Delta lines: +2 -2 =================================================================== --- TOVSolverHot/param.ccl 2010-12-04 21:45:57 UTC (rev 24) +++ TOVSolverHot/param.ccl 2010-12-04 21:53:41 UTC (rev 25) @@ -24,12 +24,12 @@ 0.0:* :: "Must be positive" } 0.0

-CCTK_REAL TOV_Gamma[10] "The polytropic constant in P = K rho^Gamma" +CCTK_REAL TOV_Gamma "The polytropic constant in P = K rho^Gamma" { 1.0: :: "The physical range at high Lorentz factors is [1,2], but otherwise higher values of gamma can also be used" } 2.0

-CCTK_REAL TOV_K[10] "The polytropic constant in P = K rho^Gamma" +CCTK_REAL TOV_K "The polytropic constant in P = K rho^Gamma" { (0.0:* :: "Greater than 0" } 100.0

Directory: /TOVSolverHot/src/ =============================

File [modified]: external.inc Delta lines: +36 -36 =================================================================== --- TOVSolverHot/src/external.inc 2010-12-04 21:45:57 UTC (rev 24) +++ TOVSolverHot/src/external.inc 2010-12-04 21:53:41 UTC (rev 25) @@ -87,8 +87,8 @@ &r_to_star, TOV_Surface[star], &press, &phi, &r); press = (press > 0.0) ? press : 0.0; - rho = pow(press/TOV_K[star], 1.0/TOV_Gamma[star]); - eps = (rho > 0.0) ? press/(TOV_Gamma[star] - 1.0) / rho : 0.0; + rho = pow(press/TOV_K, 1.0/TOV_Gamma); + eps = (rho > 0.0) ? press/(TOV_Gamma - 1.0) / rho : 0.0; gxx = r / (r_to_star + 1.0e-30) * r / (r_to_star + 1.0e-30); /* velocity components as in gr-qc/9811015 */ w_lorentz_2 = 1. / (1. - gxx*TOV_Velocity_x[star]*TOV_Velocity_x[star] @@ -97,7 +97,7 @@ if (rho > 0.0) /* source[i]=gxx*gxx*rho*(1.0+eps); */ source[i]=gxx*gxx* (w_lorentz_2*rho + - (w_lorentz_2*TOV_Gamma[star]/(TOV_Gamma[star]-1.)-1.)* + (w_lorentz_2*TOV_Gamma/(TOV_Gamma-1.)-1.)* press); else source[i]=0.0; @@ -156,14 +156,14 @@ &(TOV_rbar_1d[star_i]), &(TOV_r_1d[star_i]), &r_to_star, TOV_Surface[star], &press, &phi, &r); - rho = pow(press/TOV_K[star], 1.0/TOV_Gamma[star]); + rho = pow(press/TOV_K, 1.0/TOV_Gamma); gxx = r / (r_to_star + 1.0e-30) * r / (r_to_star + 1.0e-30); psip = pow(gxx, TOV_Momentum_Psi_Power/4.); /* velocity components as in gr-qc/9811015 */ w_lorentz_2 = 1./(1. - gxx*TOV_Velocity_x[star]*TOV_Velocity_x[star] - gxx*TOV_Velocity_y[star]*TOV_Velocity_y[star] - gxx*TOV_Velocity_z[star]*TOV_Velocity_z[star]); - rho_eps = w_lorentz_2 * (rho + TOV_Gamma[star]/(TOV_Gamma[star]-1.) * press); + rho_eps = w_lorentz_2 * (rho + TOV_Gamma/(TOV_Gamma-1.) * press); switch(dir) { case 0: source[i]=psip*rho_eps*gxx*TOV_Velocity_x[star]; break; @@ -228,9 +228,9 @@ &(TOV_rbar_1d[star_i]), &(TOV_r_1d[star_i]), &r_to_star, TOV_Surface[star], &press, &phi, &r); - rho = pow(press/TOV_K[star], 1.0/TOV_Gamma[star]); + rho = pow(press/TOV_K, 1.0/TOV_Gamma); if (rho>0.0) - eps = press/(TOV_Gamma[star] - 1.0) / rho; + eps = press/(TOV_Gamma - 1.0) / rho; else eps = 0.0; gxx = r / (r_to_star + 1.0e-30) * r / (r_to_star + 1.0e-30); @@ -460,11 +460,11 @@ CCTK_REAL f, df; count++; vx = mom_source[0][i3D]/my_psi4/ - (source[i3D]/my_psi4/my_psi4+TOV_K[0] * pow(rhonew, TOV_Gamma[0])); + (source[i3D]/my_psi4/my_psi4+TOV_K * pow(rhonew, TOV_Gamma[0])); vy = mom_source[1][i3D]/my_psi4/ - (source[i3D]/my_psi4/my_psi4+TOV_K[0] * pow(rhonew, TOV_Gamma[0])); + (source[i3D]/my_psi4/my_psi4+TOV_K * pow(rhonew, TOV_Gamma[0])); vz = mom_source[2][i3D]/my_psi4/ - (source[i3D]/my_psi4/my_psi4+TOV_K[0] * pow(rhonew, TOV_Gamma[0])); + (source[i3D]/my_psi4/my_psi4+TOV_K * pow(rhonew, TOV_Gamma[0])); v_2 = my_psi4*(vx*vx+vy*vy+vz*vz); if (count > 100) CCTK_VWarn(count<110, __LINE__, __FILE__, CCTK_THORNSTRING, \ @@ -475,36 +475,36 @@ /* velocity components as in gr-qc/9811015 */ w_lorentz_2 = 1./(1.- v_2); /* - f = TOV_K[0] / (TOV_Gamma[0]-1.0) * pow(rhonew, TOV_Gamma[0]) + + f = TOV_K / (TOV_Gamma-1.0) * pow(rhonew, TOV_Gamma) + rhonew - source[i3D]; - df = TOV_K[0] * TOV_Gamma[0] / (TOV_Gamma[0]-1.0) * - pow(rhonew, TOV_Gamma[0]-1.0) + 1.0; + df = TOV_K * TOV_Gamma / (TOV_Gamma-1.0) * + pow(rhonew, TOV_Gamma-1.0) + 1.0; */ f = my_psi4*my_psi4*( w_lorentz_2 * rhonew + - (w_lorentz_2*TOV_Gamma[0]/(TOV_Gamma[0]-1.)-1.) * - TOV_K[0] * pow(rhonew, TOV_Gamma[0])) - source[i3D]; + (w_lorentz_2*TOV_Gamma/(TOV_Gamma-1.)-1.) * + TOV_K * pow(rhonew, TOV_Gamma)) - source[i3D]; /* df= w_lorentz_2 + - (w_lorentz_2*TOV_Gamma[0]/(TOV_Gamma[0]-1.)-1.) * - TOV_K[0] * TOV_Gamma[0] * pow(rhonew, TOV_Gamma[0]-1.); + (w_lorentz_2*TOV_Gamma/(TOV_Gamma-1.)-1.) * + TOV_K * TOV_Gamma * pow(rhonew, TOV_Gamma-1.); */ /* d_w_lorentz_2/drhonew */ CCTK_REAL d_w_lorentz_2 = -2 * w_lorentz_2*w_lorentz_2 * - v_2 * TOV_K[0]*TOV_Gamma[0]*pow(rhonew, TOV_Gamma[0]-1) / - (source[i3D]/my_psi4/my_psi4 + TOV_K[0]*pow(rhonew, TOV_Gamma[0])); + v_2 * TOV_K*TOV_Gamma*pow(rhonew, TOV_Gamma-1) / + (source[i3D]/my_psi4/my_psi4 + TOV_K*pow(rhonew, TOV_Gamma)); df = my_psi4*my_psi4*( d_w_lorentz_2*rhonew + w_lorentz_2 + - d_w_lorentz_2*TOV_Gamma[0]/(TOV_Gamma[0]-1.)* - TOV_K[0] * pow(rhonew, TOV_Gamma[0]) + - (w_lorentz_2*TOV_Gamma[0]/(TOV_Gamma[0]-1.)-1.) * - TOV_K[0] * TOV_Gamma[0] * pow(rhonew, TOV_Gamma[0]-1.)); + d_w_lorentz_2*TOV_Gamma/(TOV_Gamma-1.)* + TOV_K * pow(rhonew, TOV_Gamma) + + (w_lorentz_2*TOV_Gamma/(TOV_Gamma-1.)-1.) * + TOV_K * TOV_Gamma * pow(rhonew, TOV_Gamma-1.)); /* df= 1. + my_psi4*v_2/(rhonew*rhonew) + - ((TOV_Gamma[0]+(2.-TOV_Gamma[0])*my_psi4*v_2/(rhonew*rhonew)) - / (TOV_Gamma[0]-1.) - 1.) - * TOV_K[0] * TOV_Gamma[0] * pow(rhonew, TOV_Gamma[0]-1.); + ((TOV_Gamma+(2.-TOV_Gamma)*my_psi4*v_2/(rhonew*rhonew)) + / (TOV_Gamma-1.) - 1.) + * TOV_K * TOV_Gamma * pow(rhonew, TOV_Gamma-1.); */

rhoold = rhonew; @@ -516,22 +516,22 @@ if (count==1) printf("Fast rescale! %g %g %g %g\n", my_psi4, my_psi4, mom_source[0][i3D]/my_psi4/ - (source[i3D]/my_psi4/my_psi4+TOV_K[0]* - pow(rhonew, TOV_Gamma[0])), + (source[i3D]/my_psi4/my_psi4+TOV_K* + pow(rhonew, TOV_Gamma)), mom_source[0][i3D]/my_psi4/ - (source[i3D]/my_psi4/my_psi4+TOV_K[0]* - pow(rhonew, TOV_Gamma[0]))); + (source[i3D]/my_psi4/my_psi4+TOV_K* + pow(rhonew, TOV_Gamma))); #else newton_raphson(&vx,&vy,&vz,&rhoold,&rhonew,&w_lorentz_2,&v_2, - TOV_Gamma[0], TOV_K[0], + TOV_Gamma, TOV_K, source[i3D], mom_source[0][i3D],mom_source[1][i3D],mom_source[2][i3D], my_psi4, x[i3D], y[i3D], z[i3D], rho[i3D]); #endif rho[i3D] = rhonew; - press[i3D] = TOV_K[0] * pow(rhonew, TOV_Gamma[0]); - eps[i3D] = TOV_K[0] * pow(rhonew, TOV_Gamma[0] - 1.0) / - (TOV_Gamma[0] - 1.0); + press[i3D] = TOV_K * pow(rhonew, TOV_Gamma); + eps[i3D] = TOV_K * pow(rhonew, TOV_Gamma - 1.0) / + (TOV_Gamma - 1.0); sqrt_det = pow(my_psi4, 1.5); w_lorentz[i3D] = sqrt(w_lorentz_2); dens[i3D] = sqrt_det * w_lorentz[i3D] * rho[i3D]; @@ -558,11 +558,11 @@ "%.15g %.15g %.15g %.15g\n", /*7*/ gxx[i3D]*gxx[i3D]*( /* ham - source*/ w_lorentz_2*rho[i3D]+ - (w_lorentz_2*TOV_Gamma[0]/(TOV_Gamma[0]-1.)-1.)*press[i3D]), + (w_lorentz_2*TOV_Gamma/(TOV_Gamma-1.)-1.)*press[i3D]), gxx[i3D]*gxx[i3D]*( /* ham - source */ w_lorentz_2*(rho[i3D]*(1.+eps[i3D])+press[i3D])-press[i3D]), (w_lorentz_2 * rho[i3D] + /* momx - source */ - w_lorentz_2*TOV_Gamma[0]/(TOV_Gamma[0]-1.) * press[i3D])* + w_lorentz_2*TOV_Gamma/(TOV_Gamma-1.) * press[i3D])* gxx[i3D]*vel0[i3D], w_lorentz_2 * (rho[i3D]* /* momx - source */ (1.+eps[i3D])+ press[i3D])*

File [modified]: tov.c Delta lines: +17 -17 =================================================================== --- TOVSolverHot/src/tov.c 2010-12-04 21:45:57 UTC (rev 24) +++ TOVSolverHot/src/tov.c 2010-12-04 21:53:41 UTC (rev 25) @@ -262,8 +262,8 @@ TOV_C_fill(&(TOV_rprop_1d[star_i]), TOV_Num_Radial, 0.0);

/* set start values */ - TOV_press_1d[star_i] = TOV_K[star] * - pow(rho_central, TOV_Gamma[star]); + TOV_press_1d[star_i] = TOV_K * + pow(rho_central, TOV_Gamma); /* TOV_r_1d [star_i] = LOCAL_TINY; TOV_rbar_1d [star_i] = LOCAL_TINY;*/

@@ -295,25 +295,25 @@ TOV_C_fill(source_data, 6, 0.0);

TOV_C_Source_RHS(TOV_r_1d[i], - TOV_K[star], TOV_Gamma[star], + TOV_K, TOV_Gamma, in_data, source_data);

RKLOOP k1[rk] = TOV_dr[star] * source_data[rk]; RKLOOP in_data[rk] = old_data[rk] + 0.5 * k1[rk]; TOV_C_Source_RHS(TOV_r_1d[i]+ 0.5 * TOV_dr[star], - TOV_K[star], TOV_Gamma[star], + TOV_K, TOV_Gamma, in_data, source_data);

RKLOOP k2[rk] = TOV_dr[star] * source_data[rk]; RKLOOP in_data[rk] = old_data[rk] + 0.5 * k2[rk]; TOV_C_Source_RHS(TOV_r_1d[i]+ 0.5 * TOV_dr[star], - TOV_K[star], TOV_Gamma[star], + TOV_K, TOV_Gamma, in_data, source_data);

RKLOOP k3[rk] = TOV_dr[star] * source_data[rk]; RKLOOP in_data[rk] = old_data[rk] + k3[rk]; TOV_C_Source_RHS(TOV_r_1d[i]+ TOV_dr[star], - TOV_K[star], TOV_Gamma[star], + TOV_K, TOV_Gamma, in_data, source_data); RKLOOP k4[rk] = TOV_dr[star] * source_data[rk]; RKLOOP new_data[rk] = old_data[rk] + (k1[rk] + k4[rk] + 2.0 * (k2[rk] + k3[rk])) /6.0; @@ -329,7 +329,7 @@ if (TOV_press_1d[i+1] < 0.0) TOV_press_1d[i+1] = 0.0;

- local_rho = pow(TOV_press_1d[i+1] / TOV_K[star], 1.0 / TOV_Gamma[star]); + local_rho = pow(TOV_press_1d[i+1] / TOV_K, 1.0 / TOV_Gamma);

/* scan for the surface */ if ( (local_rho <= 0.0) || @@ -375,7 +375,7 @@ CCTK_VInfo(CCTK_THORNSTRING, "Information about the TOVs used:"); CCTK_VInfo("", "TOV radius mass bary_mass mass(g) cent.rho rho(cgi) K K(cgi) Gamma"); for (i=0; i<TOV_Num_TOVs; i++) - if (fabs(TOV_Gamma[i] - 2.0) < LOCAL_TINY) + if (fabs(TOV_Gamma - 2.0) < LOCAL_TINY) CCTK_VInfo(""," %d %8g %8g %8g %8.3g %8g %8.3g %8g %8.3g %8g", i+1, TOV_R_Surface[i], TOV_m_1d[(i+1)*TOV_Num_Radial-1], @@ -385,11 +385,11 @@ TOV_Rho_Central[i]/pow(CONSTANT_G_cgi,3.0)/ pow(CONSTANT_Msolar_cgi,2.0)* pow(CONSTANT_c_cgi,6.0), - TOV_K[i], - TOV_K[i]*pow(CONSTANT_G_cgi,3.0)* + TOV_K, + TOV_K*pow(CONSTANT_G_cgi,3.0)* pow(CONSTANT_Msolar_cgi,2.0)/ pow(CONSTANT_c_cgi,4.0), - TOV_Gamma[i]); + TOV_Gamma); else CCTK_VInfo(""," %d %8g %8g %8.3g %8g %8.3g %8g %8g", i+1, TOV_R_Surface[i], @@ -399,7 +399,7 @@ TOV_Rho_Central[i]/pow(CONSTANT_G_cgi,3.0)/ pow(CONSTANT_Msolar_cgi,2.0)* pow(CONSTANT_c_cgi,6.0), - TOV_K[i], TOV_Gamma[i]); + TOV_K, TOV_Gamma);

}

@@ -643,17 +643,17 @@

/* is some perturbation wanted? */ if (Perturb[star] == 0) - rho_point[star] = pow(press_point[star]/TOV_K[star], - 1.0/TOV_Gamma[star]); + rho_point[star] = pow(press_point[star]/TOV_K, + 1.0/TOV_Gamma); else - rho_point[star] = pow(press_point[star]/TOV_K[star], - 1.0/TOV_Gamma[star]) * + rho_point[star] = pow(press_point[star]/TOV_K, + 1.0/TOV_Gamma) * (1.0 + Pert_Amplitude[star] * cos(PI/2.0 * r[i3D] / TOV_R_Surface[star]));

if (rho_point[star] > local_tiny) - eps_point[star] = press_point[star] / (TOV_Gamma[star] - 1.0) + eps_point[star] = press_point[star] / (TOV_Gamma - 1.0) / rho_point[star]; else eps_point[star] = 0.0;