User: cott Date: 2010/12/04 06:09 PM
Modified: /TOVSolverHot/ param.ccl /TOVSolverHot/src/ tov.c
Log: * this TOV solver can not make microphysical NSs
File Changes:
Directory: /TOVSolverHot/ =========================
File [modified]: param.ccl Delta lines: +15 -0 =================================================================== --- TOVSolverHot/param.ccl 2010-12-04 22:52:47 UTC (rev 26) +++ TOVSolverHot/param.ccl 2010-12-05 00:09:30 UTC (rev 27) @@ -155,11 +155,23 @@ 2 :: "entropy" } 1
+CCTK_REAL nuc_eos_rho_min "Minimum density to use with nuc_eos [Cactus]" +{ + 0.0:* :: "Better be positive or zero" +} 1.0e-13 + CCTK_REAL nuc_eos_temperature "Temperature to use when using nuc_eos [MeV]" { 0.1:* :: "Below 0.1 MeV tables get bad" } 1.0e0
+ +CCTK_REAL nuc_eos_ye "Y_e to use when using nuc_eos" +{ + 0.1:* :: "Below 0.1 things get bad" +} 0.1e0 + + CCTK_REAL nuc_eos_entropy "Entropy to use when using nuc_eos [k_B/baryon]" { 0.3:* :: "Below 0.3 k_B/baryon things get dicey" @@ -179,6 +191,9 @@ "tov" :: "TOV star initial hydrobase variables" }
+USES KEYWORD Y_e_evolution_method +USES KEYWORD temperature_evolution_method + shares:admbase
EXTENDS KEYWORD initial_data
Directory: /TOVSolverHot/src/ =============================
File [modified]: tov.c Delta lines: +164 -29 =================================================================== --- TOVSolverHot/src/tov.c 2010-12-04 22:52:47 UTC (rev 26) +++ TOVSolverHot/src/tov.c 2010-12-05 00:09:30 UTC (rev 27) @@ -159,10 +159,18 @@ CCTK_REAL xye; CCTK_REAL xent; if(eoskey==4) { - CCTK_WARN(0,"fix me"); + xtemp = nuc_eos_temperature; + xent = nuc_eos_entropy; + xye = nuc_eos_ye; + if(nuc_eos_keytemp==2) { + CCTK_WARN(0,"nuc_eos_keytemp=2 not yet handled"); + } else { + keytemp = 1; + } } else { keytemp = 0; } + rho = TOV_Rho_Central[0]; EOS_Omni_RhoFromPressEpsTempEnt(eoskey, keytemp, TOVSolverHot_eos_rf_prec, @@ -176,6 +184,15 @@ &keyerr, &anyerr);
+ + if(anyerr) { + CCTK_VInfo(CCTK_THORNSTRING, "EOS problem in TOV 1: %d %d", + keyerr,anyerr); + CCTK_VInfo(CCTK_THORNSTRING, "%15.6E %15.6E %15.6E %15.6E %15.6E", + rho,eps,xtemp,xye,press); + CCTK_WARN(0,"aborting"); + } + } mu = rho * (1.0 + eps);
@@ -228,6 +245,7 @@ k1[NUMVARS], k2[NUMVARS], k3[NUMVARS], k4[NUMVARS]; CCTK_REAL Surface_Mass, factor, local_rho; CCTK_REAL rho_central; + CCTK_INT have_temp_ye;
LOCAL_TINY = 1.0e-20;
@@ -244,6 +262,15 @@ assert(TOV_mbary_1d!=0); assert(TOV_rprop_1d!=0);
+ if(use_EOS_Omni && eoskey == 4) { + if(CCTK_EQUALS(Y_e_evolution_method,"none") || + CCTK_EQUALS(temperature_evolution_method,"none")) { + CCTK_WARN(0,"temperature and/or Y_e evolution not active"); + } else { + have_temp_ye = 1; + } + } + /* do it for all stars */ for (star=0; star < TOV_Num_TOVs; star++) { @@ -305,7 +332,14 @@ CCTK_INT n = 1; CCTK_REAL xeps,xtemp,xye,xent; if(eoskey==4) { - CCTK_WARN(0,"fix me"); + xtemp = nuc_eos_temperature; + xent = nuc_eos_entropy; + xye = nuc_eos_ye; + if(nuc_eos_keytemp==2) { + CCTK_WARN(0,"nuc_eos_keytemp=2 not yet handled"); + } else { + keytemp = 1; + } } else { keytemp = 0; } @@ -320,6 +354,14 @@ &(TOV_press_1d[star_i]), &keyerr, &anyerr); + if(anyerr) { + CCTK_VInfo(CCTK_THORNSTRING, "EOS problem in TOV 2: %d", + keyerr); + CCTK_VInfo(CCTK_THORNSTRING, "%15.6E %15.6E %15.6E %15.6E %15.6E", + rho_central,eps,xtemp,xye,TOV_press_1d[star_i]); + CCTK_WARN(0,"aborting"); + } + }
@@ -398,10 +440,18 @@ CCTK_INT n = 1; CCTK_REAL xeps,xtemp,xye,xent; if(eoskey==4) { - CCTK_WARN(0,"fix me"); + xye = nuc_eos_ye; + xtemp = nuc_eos_temperature; + xent = nuc_eos_entropy; + if(nuc_eos_keytemp==2) { + CCTK_WARN(0,"nuc_eos_keytemp=2 not yet handled"); + } else { + keytemp = 1; + } } else { keytemp = 0; } + local_rho = TOV_Rho_Central[0]; EOS_Omni_RhoFromPressEpsTempEnt(eoskey, keytemp, TOVSolverHot_eos_rf_prec, @@ -415,16 +465,36 @@ &keyerr, &anyerr);
+ if(anyerr) { + CCTK_VInfo(CCTK_THORNSTRING, "EOS problem in TOV 3: %d", + keyerr); + CCTK_VInfo(CCTK_THORNSTRING, "%15.6E %15.6E %15.6E %15.6E %15.6E", + local_rho,xeps,xtemp,xye,TOV_press_1d[i+1]); + CCTK_WARN(0,"aborting"); + } + // CCTK_VInfo(CCTK_THORNSTRING, "%6d %15.6E %15.6E", + // i,local_rho,TOV_rbar_1d[i]); + }
/* scan for the surface */ - if ( (local_rho <= 0.0) || - (TOV_press_1d[i+1] <= 0.0) ) - { - TOV_Surface[star] = TOV_r_1d[i]; - TOV_R_Surface[star] = TOV_rbar_1d[i]; - TOV_RProp_Surface[star] = TOV_rprop_1d[i]; - TOV_Surface_Index = i; + if(!(use_EOS_Omni && eoskey==4)) { + if ( (local_rho <= 0.0) || + (TOV_press_1d[i+1] <= 0.0) ) + { + TOV_Surface[star] = TOV_r_1d[i]; + TOV_R_Surface[star] = TOV_rbar_1d[i]; + TOV_RProp_Surface[star] = TOV_rprop_1d[i]; + TOV_Surface_Index = i; + } + } else { + if ( (local_rho <= nuc_eos_rho_min) ) + { + TOV_Surface[star] = TOV_r_1d[i]; + TOV_R_Surface[star] = TOV_rbar_1d[i]; + TOV_RProp_Surface[star] = TOV_rprop_1d[i]; + TOV_Surface_Index = i; + } } } if (TOV_Surface[star] < 0.0) @@ -476,16 +546,31 @@ pow(CONSTANT_c_cgi,4.0), TOV_Gamma); } else { - CCTK_VInfo("", "TOV radius mass bary_mass cent.rho rho(cgs)"); - for (i=0; i<TOV_Num_TOVs; i++) - CCTK_VInfo(""," %d %8g %8g %8g %8.3g %8g", - i+1, TOV_R_Surface[i], - TOV_m_1d[(i+1)*TOV_Num_Radial-1], - TOV_mbary_1d[(i+1)*TOV_Num_Radial-1], - TOV_Rho_Central[i], - TOV_Rho_Central[i]/pow(CONSTANT_G_cgi,3.0)/ - pow(CONSTANT_Msolar_cgi,2.0)* + if(eoskey != 4 ) { + CCTK_VInfo("", "TOV radius mass bary_mass cent.rho rho(cgs)"); + for (i=0; i<TOV_Num_TOVs; i++) + CCTK_VInfo(""," %d %8g %8g %8g %8.3g %8g", + i+1, TOV_R_Surface[i], + TOV_m_1d[(i+1)*TOV_Num_Radial-1], + TOV_mbary_1d[(i+1)*TOV_Num_Radial-1], + TOV_Rho_Central[i], + TOV_Rho_Central[i]/pow(CONSTANT_G_cgi,3.0)/ + pow(CONSTANT_Msolar_cgi,2.0)* pow(CONSTANT_c_cgi,6.0)); + } else { + CCTK_VInfo(CCTK_THORNSTRING, "Made a hot polytrope: T = %5.2f MeV", + nuc_eos_temperature); + CCTK_VInfo(CCTK_THORNSTRING, "TOV radius mass bary_mass cent.rho rho(cgs)"); + for (i=0; i<TOV_Num_TOVs; i++) + CCTK_VInfo(CCTK_THORNSTRING," %d %8g %8g %8g %8.3g %8g", + i+1, TOV_R_Surface[i], + TOV_m_1d[(i+1)*TOV_Num_Radial-1], + TOV_mbary_1d[(i+1)*TOV_Num_Radial-1], + TOV_Rho_Central[i], + TOV_Rho_Central[i]/pow(CONSTANT_G_cgi,3.0)/ + pow(CONSTANT_Msolar_cgi,2.0)* + pow(CONSTANT_c_cgi,6.0)); + } }
} @@ -593,7 +678,8 @@ DECLARE_CCTK_PARAMETERS
CCTK_REAL *press_point, *rho_point, *eps_point, - *mu_point, *phi_point, *r_point; + *mu_point, *phi_point, *r_point, *temperature_point, + *ye_point,*entropy_point;
CCTK_INT LSH_MAX_I;
@@ -625,14 +711,18 @@ assert(TOV_m_1d!=0);
/* allocate local arrays */ - r_to_star = (CCTK_REAL *) calloc (TOV_Num_TOVs, sizeof(CCTK_REAL)); - press_point = (CCTK_REAL *) calloc (TOV_Num_TOVs, sizeof(CCTK_REAL)); - rho_point = (CCTK_REAL *) calloc (TOV_Num_TOVs, sizeof(CCTK_REAL)); - eps_point = (CCTK_REAL *) calloc (TOV_Num_TOVs, sizeof(CCTK_REAL)); - mu_point = (CCTK_REAL *) calloc (TOV_Num_TOVs, sizeof(CCTK_REAL)); - phi_point = (CCTK_REAL *) calloc (TOV_Num_TOVs, sizeof(CCTK_REAL)); - r_point = (CCTK_REAL *) calloc (TOV_Num_TOVs, sizeof(CCTK_REAL)); + r_to_star = (CCTK_REAL *) calloc (TOV_Num_TOVs, sizeof(CCTK_REAL)); + press_point = (CCTK_REAL *) calloc (TOV_Num_TOVs, sizeof(CCTK_REAL)); + rho_point = (CCTK_REAL *) calloc (TOV_Num_TOVs, sizeof(CCTK_REAL)); + eps_point = (CCTK_REAL *) calloc (TOV_Num_TOVs, sizeof(CCTK_REAL)); + temperature_point = (CCTK_REAL *) calloc (TOV_Num_TOVs, sizeof(CCTK_REAL)); + ye_point = (CCTK_REAL *) calloc (TOV_Num_TOVs, sizeof(CCTK_REAL)); + entropy_point = (CCTK_REAL *) calloc (TOV_Num_TOVs, sizeof(CCTK_REAL)); + mu_point = (CCTK_REAL *) calloc (TOV_Num_TOVs, sizeof(CCTK_REAL)); + phi_point = (CCTK_REAL *) calloc (TOV_Num_TOVs, sizeof(CCTK_REAL)); + r_point = (CCTK_REAL *) calloc (TOV_Num_TOVs, sizeof(CCTK_REAL));
+ /* clear initial data */ if (TOV_Clear_Initial_Data > 0 && !(TOV_Use_Old_Initial_Data)) { @@ -692,6 +782,12 @@ TOV_C_fill(velx, LSH_MAX_I+1, 0.0); TOV_C_fill(vely, LSH_MAX_I+1, 0.0); TOV_C_fill(velz, LSH_MAX_I+1, 0.0); + if(eoskey == 4) { + TOV_C_fill(temperature, LSH_MAX_I+1, 0.0); + TOV_C_fill(entropy, LSH_MAX_I+1, 0.0); + TOV_C_fill(Y_e, LSH_MAX_I+1, 0.0); + } + } /* use the fast interpolation? only useful for testing this */ if (TOV_Fast_Interpolation == 0) @@ -741,7 +837,14 @@ CCTK_INT n = 1; CCTK_REAL xeps,xtemp,xye,xent; if(eoskey==4) { - CCTK_WARN(0,"fix me"); + xtemp = nuc_eos_temperature; + xent = nuc_eos_entropy; + xye = nuc_eos_ye; + if(nuc_eos_keytemp==2) { + CCTK_WARN(0,"nuc_eos_keytemp=2 not yet handled"); + } else { + keytemp = 1; + } } else { keytemp = 0; } @@ -758,6 +861,18 @@ &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"); + } + + temperature_point[star] = xtemp; + entropy_point[star] = xent; + ye_point[star] = xye; + } else if(!use_EOS_Omni) { @@ -899,7 +1014,12 @@ rho[i3D] = TOV_Atmosphere[star]; eps[i3D] = eps_point[star]; press[i3D] = press_point[star]; - + + if(eoskey==4) { + temperature[i3D] = temperature_point[star]; + entropy[i3D] = entropy_point[star]; + Y_e[i3D] = ye_point[star]; + } } else if (CCTK_EQUALS(TOV_Combine_Method, "average")) { @@ -935,6 +1055,12 @@ rho[i3D] += rho_point[star_i]; eps[i3D] += eps_point[star_i]; press[i3D] += press_point[star_i]; + if(eoskey==4) { + temperature[i3D] += temperature_point[star]; + entropy[i3D] += entropy_point[star]; + Y_e[i3D] += ye_point[star]; + } + /* we still have to know if we are inside one star - and which */ if ( (rho_point[star_i] > max_rho) && (rho_point[star_i] > TOV_Atmosphere[star_i])) @@ -1005,6 +1131,11 @@ (1.0 + eps[i3D]) + press[i3D]*(w_lorentz[i3D]*w_lorentz[i3D]-1.0) ) - dens[i3D]; + + if(eoskey==4) { + Y_e_con[i3D] = Y_e[i3D] * dens[i3D]; + } + } /* if used, recalculate the derivatives of the conformal factor */ if (*conformal_state > 1) @@ -1095,6 +1226,10 @@ free(mu_point); free(phi_point); free(r_point); + free(temperature_point); + free(entropy_point); + free(ye_point); + }
void TOV_Prepare_Fake_Evolution(CCTK_ARGUMENTS)