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

4 Dec 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 */

    