[Commits] [svn:einsteintoolkit] incoming/TOVSolverHot/ (Rev. 26)
cott at tapir.caltech.edu
cott at tapir.caltech.edu
Sat Dec 4 16:52:47 CST 2010
User: cott
Date: 2010/12/04 04:52 PM
Modified:
/TOVSolverHot/
interface.ccl, param.ccl
/TOVSolverHot/src/
external.inc, tov.c
Log:
* this now works for polytropic NSs using EOS omni.
File Changes:
Directory: /TOVSolverHot/
=========================
File [modified]: interface.ccl
Delta lines: +14 -12
===================================================================
--- TOVSolverHot/interface.ccl 2010-12-04 21:53:41 UTC (rev 25)
+++ TOVSolverHot/interface.ccl 2010-12-04 22:52:47 UTC (rev 26)
@@ -91,17 +91,19 @@
USES FUNCTION EOS_Omni_press
-void FUNCTION EOS_Omni_RhoFromPress(CCTK_INT IN eoskey, \
- CCTK_INT IN havetemp, \
- CCTK_REAL IN rf_precision, \
- CCTK_INT IN npoints, \
- CCTK_REAL INOUT ARRAY rho, \
- CCTK_REAL INOUT ARRAY eps, \
- CCTK_REAL IN ARRAY temp, \
- CCTK_REAL IN ARRAY ye, \
- CCTK_REAL IN ARRAY press, \
- CCTK_REAL OUT ARRAY xeps, \
- CCTK_INT OUT ARRAY keyerr, \
+
+void FUNCTION EOS_Omni_RhoFromPressEpsTempEnt( \
+ CCTK_INT IN eoskey, \
+ CCTK_INT IN havetemp, \
+ CCTK_REAL IN rf_precision, \
+ CCTK_INT IN npoints, \
+ CCTK_REAL OUT ARRAY rho, \
+ CCTK_REAL INOUT ARRAY eps, \
+ CCTK_REAL INOUT ARRAY temp, \
+ CCTK_REAL INOUT ARRAY entropy, \
+ CCTK_REAL IN ARRAY ye, \
+ CCTK_REAL IN ARRAY press, \
+ CCTK_INT OUT ARRAY keyerr, \
CCTK_INT OUT anyerr)
-USES FUNCTION EOS_Omni_RhoFromPress
+USES FUNCTION EOS_Omni_RhoFromPressEpsTempEnt
File [modified]: param.ccl
Delta lines: +41 -0
===================================================================
--- TOVSolverHot/param.ccl 2010-12-04 21:53:41 UTC (rev 25)
+++ TOVSolverHot/param.ccl 2010-12-04 22:52:47 UTC (rev 26)
@@ -131,6 +131,47 @@
".*" :: "Any filename, not used if empty"
} ""
+########################################################
+# EOS Omni Parameters
+########################################################
+BOOLEAN use_EOS_Omni "should we use omni?"
+{
+} no
+
+CCTK_REAL TOVSolverHot_eos_rf_prec "precision for root finder"
+{
+ 1.0e-12:* :: "anything above crazy"
+} 1.0e-10
+
+CCTK_INT eoskey "Which EOS are we using"
+{
+ 1 :: "Polytropic"
+ 4 :: "nuc_eos"
+} 1
+
+CCTK_INT nuc_eos_keytemp "When using nuc_eos: fix temperature or entropy?"
+{
+ 1 :: "temperature"
+ 2 :: "entropy"
+} 1
+
+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_entropy "Entropy to use when using nuc_eos [k_B/baryon]"
+{
+ 0.3:* :: "Below 0.3 k_B/baryon things get dicey"
+} 1.0e0
+
+
+
+
+########################################################
+# Other Parameters
+########################################################
+
shares:HydroBase
EXTENDS KEYWORD initial_hydro ""
Directory: /TOVSolverHot/src/
=============================
File [modified]: external.inc
Delta lines: +21 -0
===================================================================
--- TOVSolverHot/src/external.inc 2010-12-04 21:53:41 UTC (rev 25)
+++ TOVSolverHot/src/external.inc 2010-12-04 22:52:47 UTC (rev 26)
@@ -65,6 +65,9 @@
/* TOV_debug_input_points(size,x,y,z); */
/* loop over all points we got */
+ if(use_EOS_Omni)
+ CCTK_WARN(0,"TOV_Set_Rho_ADM does not work with EOS_Omni");
+
#pragma omp parallel for
for (CCTK_INT i=0; i<size; i++)
{
@@ -139,6 +142,9 @@
/* CCTK_INFO("Request for Momentum-Sources from TOVSolver"); */
/* loop over all points we got */
+ if(use_EOS_Omni)
+ CCTK_WARN(0,"TOV_Set_Momentum_Source does not work with EOS_Omni");
+
for (i=0; i<size; i++)
{
/* calculate the distance to the star */
@@ -209,6 +215,9 @@
CCTK_INFO("Using initial guess from TOVSolver");
/* TOV_debug_input_points(size,x,y,z); */
+ if(use_EOS_Omni)
+ CCTK_WARN(0,"TOV_Set_Initial_Guess_for_u does not work with EOS_Omni");
+
/* loop over all points we got */
for (i=0; i<size; i++)
{
@@ -251,9 +260,14 @@
CCTK_REAL my_psi4,
CCTK_REAL x, CCTK_REAL y, CCTK_REAL z, CCTK_REAL rho_orig)
{
+ DECLARE_CCTK_PARAMETERS;
CCTK_REAL lb, ub, f;
lb=0.0;
ub=2*rho_orig;
+
+ if(use_EOS_Omni)
+ CCTK_WARN(0,"bisection does not work with EOS_Omni");
+
do
{
*rhonew = (ub-lb)/2.;
@@ -289,8 +303,12 @@
CCTK_REAL my_psi4,
CCTK_REAL x, CCTK_REAL y, CCTK_REAL z, CCTK_REAL rho_orig)
{
+ DECLARE_CCTK_PARAMETERS
int count=0;
CCTK_REAL d_w_lorentz_2, f, df;
+ if(use_EOS_Omni)
+ CCTK_WARN(0,"newton_raphson does not work with EOS_Omni");
+
do
{
count++;
@@ -371,6 +389,9 @@
static int debug = 1;
FILE *debugfile = NULL;
+ if(use_EOS_Omni)
+ CCTK_WARN(0,"Rescaling Sources does not work with EOS_Omni");
+
CCTK_INFO("Rescaling Sources");
/* get GRHydro variables */
File [modified]: tov.c
Delta lines: +148 -27
===================================================================
--- TOVSolverHot/src/tov.c 2010-12-04 21:53:41 UTC (rev 25)
+++ TOVSolverHot/src/tov.c 2010-12-04 22:52:47 UTC (rev 26)
@@ -11,6 +11,7 @@
#include <cctk.h>
#include <cctk_Arguments.h>
+#include <cctk_Functions.h>
#include <cctk_Parameters.h>
#include "constants.h"
@@ -131,6 +132,8 @@
CCTK_REAL press, rho, eps, mu, m, phi, mbary, rprop;
CCTK_REAL log_rbar_over_r, r_minus_two_m;
+ DECLARE_CCTK_PARAMETERS;
+
LOCAL_TINY = 1.0e-35;
PI=4.0*atan(1.0);
@@ -143,8 +146,37 @@
mbary = old_data[4];
rprop = old_data[5];
- rho = pow(press / K, 1.0 / Gamma);
- eps = press / (Gamma - 1.0) / rho;
+ if(!use_EOS_Omni) {
+ rho = pow(press / K, 1.0 / Gamma);
+ eps = press / (Gamma - 1.0) / rho;
+ } else {
+ // EOS_Omni code
+ CCTK_INT keytemp;
+ CCTK_INT keyerr = 0;
+ CCTK_INT anyerr = 0;
+ CCTK_INT n = 1;
+ CCTK_REAL xtemp;
+ CCTK_REAL xye;
+ CCTK_REAL xent;
+ if(eoskey==4) {
+ CCTK_WARN(0,"fix me");
+ } else {
+ keytemp = 0;
+ }
+ EOS_Omni_RhoFromPressEpsTempEnt(eoskey,
+ keytemp,
+ TOVSolverHot_eos_rf_prec,
+ n,
+ &rho,
+ &eps,
+ &xtemp,
+ &xent,
+ &xye,
+ &press,
+ &keyerr,
+ &anyerr);
+
+ }
mu = rho * (1.0 + eps);
r_minus_two_m = r - 2.0 * m;
@@ -262,8 +294,35 @@
TOV_C_fill(&(TOV_rprop_1d[star_i]), TOV_Num_Radial, 0.0);
/* set start values */
- TOV_press_1d[star_i] = TOV_K *
- pow(rho_central, TOV_Gamma);
+ if(!use_EOS_Omni) {
+ TOV_press_1d[star_i] = TOV_K *
+ pow(rho_central, TOV_Gamma);
+ } else {
+ // EOS_Omni code
+ CCTK_INT keytemp;
+ CCTK_INT keyerr = 0;
+ CCTK_INT anyerr = 0;
+ CCTK_INT n = 1;
+ CCTK_REAL xeps,xtemp,xye,xent;
+ if(eoskey==4) {
+ CCTK_WARN(0,"fix me");
+ } else {
+ keytemp = 0;
+ }
+ EOS_Omni_press(eoskey,
+ keytemp,
+ TOVSolverHot_eos_rf_prec,
+ n,
+ &rho_central,
+ &xeps,
+ &xtemp,
+ &xye,
+ &(TOV_press_1d[star_i]),
+ &keyerr,
+ &anyerr);
+ }
+
+
/* TOV_r_1d [star_i] = LOCAL_TINY;
TOV_rbar_1d [star_i] = LOCAL_TINY;*/
@@ -329,8 +388,35 @@
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, 1.0 / TOV_Gamma);
+ if(!use_EOS_Omni) {
+ local_rho = pow(TOV_press_1d[i+1] / TOV_K, 1.0 / TOV_Gamma);
+ } else {
+ // EOS_Omni code
+ CCTK_INT keytemp;
+ CCTK_INT keyerr = 0;
+ CCTK_INT anyerr = 0;
+ CCTK_INT n = 1;
+ CCTK_REAL xeps,xtemp,xye,xent;
+ if(eoskey==4) {
+ CCTK_WARN(0,"fix me");
+ } else {
+ keytemp = 0;
+ }
+ EOS_Omni_RhoFromPressEpsTempEnt(eoskey,
+ keytemp,
+ TOVSolverHot_eos_rf_prec,
+ n,
+ &local_rho,
+ &xeps,
+ &xtemp,
+ &xent,
+ &xye,
+ &(TOV_press_1d[i+1]),
+ &keyerr,
+ &anyerr);
+ }
+
/* scan for the surface */
if ( (local_rho <= 0.0) ||
(TOV_press_1d[i+1] <= 0.0) )
@@ -373,14 +459,13 @@
CCTK_INFO("Integrated TOV equation");
/* do some info */
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 - 2.0) < LOCAL_TINY)
- CCTK_VInfo(""," %d %8g %8g %8g %8.3g %8g %8.3g %8g %8.3g %8g",
+ if(!use_EOS_Omni) {
+ CCTK_VInfo("", "TOV radius mass bary_mass cent.rho rho(cgs) K K(cgs) Gamma");
+ for (i=0; i<TOV_Num_TOVs; i++)
+ CCTK_VInfo(""," %d %8g %8g %8g %8g %8.3g %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_m_1d[(i+1)*TOV_Num_Radial-1]*CONSTANT_Msolar_cgi,
TOV_Rho_Central[i],
TOV_Rho_Central[i]/pow(CONSTANT_G_cgi,3.0)/
pow(CONSTANT_Msolar_cgi,2.0)*
@@ -390,16 +475,18 @@
pow(CONSTANT_Msolar_cgi,2.0)/
pow(CONSTANT_c_cgi,4.0),
TOV_Gamma);
- else
- CCTK_VInfo(""," %d %8g %8g %8.3g %8g %8.3g %8g %8g",
+ } 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_m_1d[(i+1)*TOV_Num_Radial-1]*CONSTANT_Msolar_cgi,
+ 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),
- TOV_K, TOV_Gamma);
+ pow(CONSTANT_c_cgi,6.0));
+ }
}
@@ -643,20 +730,54 @@
/* is some perturbation wanted? */
if (Perturb[star] == 0)
- rho_point[star] = pow(press_point[star]/TOV_K,
- 1.0/TOV_Gamma);
+ if(!use_EOS_Omni) {
+ rho_point[star] = pow(press_point[star]/TOV_K,
+ 1.0/TOV_Gamma);
+ } else {
+ // EOS_Omni code
+ CCTK_INT keytemp;
+ CCTK_INT keyerr = 0;
+ CCTK_INT anyerr = 0;
+ CCTK_INT n = 1;
+ CCTK_REAL xeps,xtemp,xye,xent;
+ if(eoskey==4) {
+ CCTK_WARN(0,"fix me");
+ } else {
+ keytemp = 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);
+
+ }
else
- 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(!use_EOS_Omni) {
+ 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]));
+ } else {
+ CCTK_WARN(0, "Perturbation handling not implemented for EOS_Omni");
+ }
- if (rho_point[star] > local_tiny)
- eps_point[star] = press_point[star] / (TOV_Gamma - 1.0)
- / rho_point[star];
- else
- eps_point[star] = 0.0;
+ if(!use_EOS_Omni) {
+ if (rho_point[star] > local_tiny)
+ eps_point[star] = press_point[star] / (TOV_Gamma - 1.0)
+ / rho_point[star];
+ else
+ eps_point[star] = 0.0;
+ } // eps already handled by EOS_Omni call above
+
mu_point[star] = rho_point[star] * (1.0 + eps_point[star]);
}
/* find out from which star we want to have the data */
More information about the Commits
mailing list