User: cott Date: 2010/12/04 08:05 PM
Modified: /TOVSolverHot/src/ tov.c
Log: * small modifications to initial guesses for rho -- this helps if one wants to make a hot NS
File Changes:
Directory: /TOVSolverHot/src/ =============================
File [modified]: tov.c Delta lines: +75 -25 =================================================================== --- TOVSolverHot/src/tov.c 2010-12-05 01:19:49 UTC (rev 29) +++ TOVSolverHot/src/tov.c 2010-12-05 02:05:49 UTC (rev 30) @@ -125,7 +125,7 @@ @history @endhistory @@*/ -void TOV_C_Source_RHS(CCTK_REAL r, CCTK_REAL K, CCTK_REAL Gamma, +void TOV_C_Source_RHS(CCTK_REAL r, CCTK_REAL K, CCTK_REAL Gamma, CCTK_REAL old_data[NUMVARS], CCTK_REAL source_data[NUMVARS]) { CCTK_REAL LOCAL_TINY, PI; @@ -170,7 +170,8 @@ } else { keytemp = 0; } - rho = TOV_Rho_Central[0]; + // guess for rho + rho = pow(press / 100.0, 1.0 / 2.0); EOS_Omni_RhoFromPressEpsTempEnt(eoskey, keytemp, TOVSolverHot_eos_rf_prec, @@ -396,7 +397,7 @@ TOV_C_fill(source_data, 6, 0.0);
TOV_C_Source_RHS(TOV_r_1d[i], - TOV_K, TOV_Gamma, + TOV_K, TOV_Gamma, in_data, source_data);
RKLOOP k1[rk] = TOV_dr[star] * source_data[rk]; @@ -835,7 +836,10 @@ CCTK_INT keyerr = 0; CCTK_INT anyerr = 0; CCTK_INT n = 1; - CCTK_REAL xeps,xtemp,xye,xent; + CCTK_REAL xeps,xtemp,xye,xent,xprs; + + + if(eoskey==4) { xtemp = nuc_eos_temperature; xent = nuc_eos_entropy; @@ -845,30 +849,76 @@ } else { keytemp = 1; } + if(press_point[star] > 0.0e0) { + // initial guess + rho_point[star] = pow(press_point[star] / 100.0, 1.0 / 2.0); + xprs = press_point[star]; + EOS_Omni_RhoFromPressEpsTempEnt(eoskey, + keytemp, + TOVSolverHot_eos_rf_prec, + n, + &(rho_point[star]), + &(eps_point[star]), + &xtemp, + &xent, + &xye, + &(press_point[star]), + &keyerr, + &anyerr); + + if(anyerr) { + CCTK_VInfo(CCTK_THORNSTRING, "EOS problem in TOV 4: %d", + keyerr); + CCTK_VInfo(CCTK_THORNSTRING, "%15.6E %15.6E %15.6E %15.6E %15.6E", + rho_point[star],eps_point[star],xtemp,xye,xprs); + CCTK_WARN(0,"aborting"); + } + } else { + rho_point[star] = nuc_eos_rho_min; + keytemp = 1; + EOS_Omni_press(eoskey, + keytemp, + TOVSolverHot_eos_rf_prec, + n, + &rho_point[star], + &eps_point[star], + &xtemp, + &xye, + &xprs, + &keyerr, + &anyerr); + if(anyerr) { + CCTK_VInfo(CCTK_THORNSTRING, "EOS problem in TOV 4a: %d", + keyerr); + CCTK_VInfo(CCTK_THORNSTRING, "%15.6E %15.6E %15.6E %15.6E %15.6E", + rho_point[star],eps_point[star],xtemp,xye,press_point[star]); + CCTK_WARN(0,"aborting"); + } + } } else { keytemp = 0; - } - // initial guess - rho_point[star] = 0.001e0; - EOS_Omni_RhoFromPressEpsTempEnt(eoskey, - keytemp, - TOVSolverHot_eos_rf_prec, - n, - &(rho_point[star]), - &(eps_point[star]), + // initial guess + rho_point[star] = pow(press_point[star] / 100.0, 1.0 / 2.0); + EOS_Omni_RhoFromPressEpsTempEnt(eoskey, + keytemp, + TOVSolverHot_eos_rf_prec, + n, + &(rho_point[star]), + &(eps_point[star]), &xtemp, - &xent, - &xye, - &(press_point[star]), - &keyerr, - &anyerr); - - if(anyerr) { - CCTK_VInfo(CCTK_THORNSTRING, "EOS problem in TOV 4: %d", - keyerr); - CCTK_VInfo(CCTK_THORNSTRING, "%15.6E %15.6E %15.6E %15.6E %15.6E", - rho_point[star],eps_point[star],xtemp,xye,press_point[star]); - CCTK_WARN(0,"aborting"); + &xent, + &xye, + &(press_point[star]), + &keyerr, + &anyerr); + + if(anyerr) { + CCTK_VInfo(CCTK_THORNSTRING, "EOS problem in TOV 4: %d", + keyerr); + CCTK_VInfo(CCTK_THORNSTRING, "%15.6E %15.6E %15.6E %15.6E %15.6E", + rho_point[star],eps_point[star],xtemp,xye,press_point[star]); + CCTK_WARN(0,"aborting"); + } }
if(rho_point[star] > nuc_eos_rho_min) {