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

cott at tapir.caltech.edu cott at tapir.caltech.edu
Sat Dec 4 18:09:30 CST 2010


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)



More information about the Commits mailing list