GRHydro: Changes to make GRHydro MHD work with nuclear/hot equation of state:

1) Switch to EOSOmni pointwise C2P routine and modify where necessary.
2) Modify Con2PrimM.F90 to allow for the evolution of temperature and adjust the wrapper routine.
3) Create EigenProblemM_hot pointwise routine and call that from HLLEM.F90 when temperature is evolved.
   Additionally adjust HLLEM where necessary.
4) Adjust InterfacesM.h to incorporate the newly created functions.
5) Fix a loop problem (not taking into account constraint transport) in PPM reconstruction of Y_e
6) Introduce Prim2ConM_hot and call this pointwise routine from Prim2ConM.F90 when temperature is evolved.
   Additionally also make this routine available to initial data routine in GRHydro_InitData
7) Adjust loops in GRHydro_PoloidalMagFieldM.F90 to not set boundary points it cannot set but instead call boundary group afterwards! Pay attention as this will not work with boundary conditions set to "none" in MHD case anymore but is the correct thing to do.
8) Allow StarMapper to extend HydroBase::initial_hydro = "starmapper".
9) Smaller fixes.

From: Philipp Moesta <pmoesta@tapir.caltech.edu>

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@410 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45

1 files changed, 35 insertions, 8 deletions

diff --git a/src/GRHydro_Con2PrimM_pt_EOSOmni.c b/src/GRHydro_Con2PrimM_pt_EOSOmni.c
index f0b6a39..4a52414 100644
--- a/src/GRHydro_Con2PrimM_pt_EOSOmni.c
+++ b/src/GRHydro_Con2PrimM_pt_EOSOmni.c
@@ -45,6 +45,7 @@
 #include <complex.h>
 
 #include "cctk.h"
+#include "cctk_Parameters.h"
 
 /* Set this to be 1 if you want debug output */
 #define DEBUG_CON2PRIMM (0)
@@ -123,7 +124,7 @@ static CCTK_REAL pressure_rho0_w(CCTK_REAL rho0, CCTK_REAL w,CCTK_REAL gammaeos)
 }
 
 void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) (
-    CCTK_INT  *handle, CCTK_INT  *keytemp, CCTK_REAL  *prec, 
+    CCTK_INT  *handle, CCTK_INT *keytemp, CCTK_REAL *prec, 
     CCTK_REAL *gamma_eos,
     CCTK_REAL *dens_in, 
     CCTK_REAL *sx_in, CCTK_REAL *sy_in, CCTK_REAL *sz_in, 
@@ -204,6 +205,8 @@ RETURN VALUE: of retval = (i*100 + j)  where
              ( used to be  5 -> failure: rho,uu <= 0   but now sets epsnegative to non-zero )
 
 **********************************************************************************/
+
+
 void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt)  ( 
     CCTK_INT  *handle, CCTK_INT *keytemp, CCTK_REAL *prec,
     CCTK_REAL *gamma_eos,
@@ -239,11 +242,13 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt)  (
   CCTK_REAL inv_sqrt_detg = 1./sqrt_detg; 
   CCTK_INT i,j, i_increase ;
 
+  DECLARE_CCTK_PARAMETERS;
+  
   struct LocGlob lg; 
   struct eosomnivars eosvars;
 
   eosvars.eoshandle = *handle;
-  printf("handle = %i\n",*handle);
+  //printf("handle = %i\n",*handle);
   eosvars.eoskeytemp = *keytemp;
   eosvars.eosprec = *prec;
   eosvars.eos_y_e[0] = *y_e_in;
@@ -271,6 +276,8 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt)  (
   fprintf(stdout," *vely        = %26.16e \n", *vely	  );    
   fprintf(stdout," *velz        = %26.16e \n", *velz	  );    
   fprintf(stdout," *epsilon     = %26.16e \n", *epsilon     );
+  fprintf(stdout," *temp_in     = %26.16e \n", *temp_in     );
+  fprintf(stdout," *y_e_in     = %26.16e \n", *y_e_in     );
   fprintf(stdout," *pressure    = %26.16e \n", *pressure    );
   fprintf(stdout," *Bx 	        = %26.16e \n", *Bx 	  );    
   fprintf(stdout," *By 	        = %26.16e \n", *By 	  );    
@@ -347,6 +354,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt)  (
   }
   if(vsq < 0. || vsq > 1. ) {
     *retval = 2.;
+    fprintf(stdout," *retval      = %26.16e \n", *retval      );
     return;
   }
 
@@ -357,7 +365,8 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt)  (
   //   i.e. you don't get positive values for dP/d(vsq) . 
   rho0 = lg.D / gamma ;
   u = (*epsilon) * rho0;
-
+  CCTK_REAL uold = u; 
+ 
   CCTK_REAL dum1,dum2;
   p = pressure_rho0_eps_eosomni(rho0,*epsilon,&dum1,&dum2,&eosvars) ;  // EOSOMNI
   //  p = pressure_rho0_u(rho0,u,gammaeos) ;  // EOS
@@ -365,6 +374,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt)  (
 
   W_last = w*gammasq ;
 
+  //fprintf(stdout," p     	 	= %26.16e \n", p	        );
 
   // Make sure that W is large enough so that v^2 < 1 : 
   i_increase = 0;
@@ -395,11 +405,13 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt)  (
   /* Problem with solver, so return denoting error before doing anything further */
   if( ((*retval) != 0.) || (W == FAIL_VAL) ) {
     *retval = *retval*100.+1.;
+    fprintf(stdout," *retval      = %26.16e \n", *retval      );
     return;
   }
   else{
     if(W <= 0. || W > W_TOO_BIG) {
       *retval = 3.;
+      fprintf(stdout," *retval      = %26.16e \n", *retval      );
       return;
     }
   }
@@ -407,6 +419,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt)  (
   // Calculate v^2:
   if( vsq >= 1. ) {
     *retval = 4.;
+    fprintf(stdout," *retval      = %26.16e \n", *retval      );
     return;
   }
 
@@ -425,13 +438,20 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt)  (
     u=x_3d[2];
     *epsilon = u/rho0;
     p = pressure_rho0_eps_eosomni(rho0,*epsilon,&dum1,&dum2,&eosvars) ;  // EOSOMNI
-    printf("%g %g %g %g\n",rho0,u,*epsilon,p);
+   // printf("%g %g %g %g\n",rho0,u,*epsilon,p);
   }
 
   // User may want to handle this case differently, e.g. do NOT return upon 
   // a negative rho/u, calculate v^i so that rho/u can be floored by other routine:
-  if( (rho0 <= 0.) || (u <= 0.) ) { 
-    *epsnegative = 1; 
+  // HOT: if( (rho0 <= 0.) || (u <= 0.) ) { 
+  if(rho0 <= 0.) { 
+//    *epsnegative = 1; 
+    fprintf(stdout," *epsnegative = %10d    \n", *epsnegative );
+    fprintf(stdout," rho0      	= %26.16e \n", rho0		); 
+    fprintf(stdout," u     	 	= %26.16e \n", u	        );
+    fprintf(stdout," W     	 	= %26.16e \n", W	        );
+    fprintf(stdout," vsq     	 	= %26.16e \n", vsq	        );
+    fprintf(stdout," uold     	 	= %26.16e \n", uold	        );
     return;
   }
 
@@ -447,6 +467,13 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt)  (
   *vely = g_o_WBsq * ( usy + QdB_o_W*(*By) ) ;
   *velz = g_o_WBsq * ( usz + QdB_o_W*(*Bz) ) ;
 
+  if (*rho <= rho_abs_min*(1.0+GRHydro_atmo_tolerance) ) {
+    *rho = rho_abs_min;
+    *velx = 0.0;
+    *vely = 0.0;
+    *velz = 0.0;
+    *w_lorentz = 1.0;
+  }
 
 #if(DEBUG_CON2PRIMM)
   fprintf(stdout,"rho          = %26.16e \n",*rho      );
@@ -696,7 +723,7 @@ static CCTK_INT threed_newton_raphson_omni( CCTK_REAL x[], struct eosomnivars *e
     x[1] += dx[1]  ; 
     x[2] += dx[2]  ; 
 
-    printf("Updating vars: %g %g %g %g %g %g\n",x[0],dx[0],x[1],dx[1],x[2],dx[2]);
+    //printf("Updating vars: %g %g %g %g %g %g\n",x[0],dx[0],x[1],dx[1],x[2],dx[2]);
 
     /****************************************/
     /* Calculate the convergence criterion */
@@ -717,7 +744,7 @@ static CCTK_INT threed_newton_raphson_omni( CCTK_REAL x[], struct eosomnivars *e
       if( x[1] > 1. ) { x[1] = dv; }
     }
 
-    if( x[2] < 0. ) {  x[2] = 0.; } 
+    // HOT: if( x[2] < 0. ) {  x[2] = 0.; } 
       
     /*****************************************************************************/
     /* If we've reached the tolerance level, then just do a few extra iterations */