[Commits] [svn:einsteintoolkit] incoming/TOVSolverHot/ (Rev. 25)
cott at tapir.caltech.edu
cott at tapir.caltech.edu
Sat Dec 4 15:53:41 CST 2010
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;
More information about the Commits
mailing list