RIT MHD development:

 Fix race conditions arising in the Con2Prim routines.
Now both GNU and Intel compilers agree on the results
for the magnetic tests.



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

3 files changed, 135 insertions, 211 deletions

diff --git a/src/GRHydro_Con2PrimM.F90 b/src/GRHydro_Con2PrimM.F90
index 7367e45..79e24eb 100644
--- a/src/GRHydro_Con2PrimM.F90
+++ b/src/GRHydro_Con2PrimM.F90
@@ -94,7 +94,7 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS)
 
   !$OMP PARALLEL DO PRIVATE(i,j,itracer,&
   !$OMP uxx, uxy, uxz, uyy, uyz, uzz, det, epsnegative, &
-  !$OMP b2,local_gam,xrho,xpress,local_K,local_pgam,sc)
+  !$OMP b2,xrho,xeps,xpress,xtemp,local_K,local_pgam,sc,keyerr,anyerr)
   do k = 1, nz 
      do j = 1, ny 
         do i = 1, nx
@@ -142,6 +142,7 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS)
               scon(i,j,k,:) = 0.d0
               vel(i,j,k,:) = 0.d0
               w_lorentz(i,j,k) = 1.d0
+              xtemp = 0.0d0
 
               call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
                    rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),keyerr,anyerr)
@@ -225,13 +226,14 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS)
               call CCTK_WARN(GRHydro_NaN_verbose+2, 'Specific internal energy just went below 0, trying polytype.')
               !$OMP END CRITICAL
               
-              xrho=10.0d0
+              xrho=1.0d0; xtemp=0.0d0; xeps=1.0d0
               call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
-                   1.d0,1.0d0,xtemp,xye,xpress,keyerr,anyerr)
-              local_K = xpress
+                   xrho,xeps,xtemp,xye,xpress,keyerr,anyerr)
+              local_K = xpress; 
 
+              xrho=10.0d0; xeps=1.0d0
               call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
-                   xrho,1.0d0,xtemp,xye,xpress,keyerr,anyerr)
+                   xrho,xeps,xtemp,xye,xpress,keyerr,anyerr)
               local_pgam=log(xpress/local_K)/log(xrho)
               sc = local_K*dens(i,j,k)
 
diff --git a/src/GRHydro_Con2PrimM_polytype_pt.c b/src/GRHydro_Con2PrimM_polytype_pt.c
index 74cd0d5..78445cd 100644
--- a/src/GRHydro_Con2PrimM_polytype_pt.c
+++ b/src/GRHydro_Con2PrimM_polytype_pt.c
@@ -71,42 +71,32 @@
 ***************************************************/
 
 /* Local Globals */
-static CCTK_REAL Bsq,QdotBsq,Qtsq,Qdotn,D,half_Bsq,rho_gm1;
-static CCTK_REAL W_for_gnr2, rho_for_gnr2, W_for_gnr2_old, rho_for_gnr2_old, drho_dW;
-static CCTK_REAL t2,t4,t7,t24,two_Bsq,t300,t400,s200,s100;
+struct LocGlob {
+  CCTK_REAL Bsq,QdotBsq,Qtsq,D;
+  CCTK_REAL W_for_gnr2, rho_for_gnr2, W_for_gnr2_old, rho_for_gnr2_old;
+  CCTK_REAL t4,t7,t24,two_Bsq,t300,t400,s200,s100;
+} ;
 
 
 // Declarations: 
-static CCTK_REAL vsq_calc(CCTK_REAL W);
-
-static CCTK_INT oned_newton_raphson( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammaeos,
+static CCTK_INT oned_newton_raphson( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammaeos, struct LocGlob *lgp,
 				     void (*funcd) (CCTK_REAL [], CCTK_REAL [], 
 						    CCTK_REAL [], CCTK_REAL [][1], 
-						    CCTK_REAL *, CCTK_REAL *, CCTK_INT, CCTK_REAL) );
+						    CCTK_REAL *, CCTK_REAL *, CCTK_INT, CCTK_REAL, 
+						    struct LocGlob *lgp) );
 
 static  void func_W( CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], CCTK_REAL [][1], 
-		     CCTK_REAL *f, CCTK_REAL *df, CCTK_INT n, CCTK_REAL gammaeos);
-static CCTK_REAL x1_of_x0(CCTK_REAL x0 ) ;
+		     CCTK_REAL *f, CCTK_REAL *df, CCTK_INT n, CCTK_REAL gammaeos, 
+		     struct LocGlob * lgp);
 
-static CCTK_INT gnr2( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammaeos,
+static CCTK_INT gnr2( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammaeos, struct LocGlob *lgp,
 			    void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], 
-					   CCTK_REAL [][1],CCTK_REAL *,CCTK_REAL *,CCTK_INT, CCTK_REAL) );
-static void func_rho(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], 
-		     CCTK_REAL jac[][1], CCTK_REAL *f, CCTK_REAL *df, CCTK_INT n, CCTK_REAL gammaeos);
-
-static CCTK_REAL eos_info(CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL *dpdw, CCTK_REAL *dpdvsq, CCTK_REAL gammaeos);
-/* pressure as a function of rho0 and u */
-static CCTK_REAL pressure_rho0_u(CCTK_REAL rho0, CCTK_REAL u, CCTK_REAL gammaeos)
-{
-  return((gammaeos - 1.)*u) ;
-}
-
-/* Pressure as a function of rho0 and w = rho0 + u + p */
-static CCTK_REAL pressure_rho0_w(CCTK_REAL rho0, CCTK_REAL w, CCTK_REAL gammaeos)
-{
-  return((gammaeos-1.)*(w - rho0)/gammaeos) ;
-}
+					   CCTK_REAL [][1],CCTK_REAL *,CCTK_REAL *,CCTK_INT, CCTK_REAL,
+					   struct LocGlob *lgp) );
 
+static void func_rho(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], 
+		     CCTK_REAL jac[][1], CCTK_REAL *f, CCTK_REAL *df, CCTK_INT n, CCTK_REAL gammaeos, 
+		     struct LocGlob *lgp);
 
 void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_Polytype_pt) (
     CCTK_INT  *handle, CCTK_REAL *gamma_eos, 
@@ -220,8 +210,12 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_Polytype_pt)  (
   CCTK_REAL detg = (*det);
   CCTK_REAL sqrt_detg = sqrt(detg);
   CCTK_REAL inv_sqrt_detg = 1./sqrt_detg; 
+  CCTK_REAL t2, rho_gm1;
   CCTK_INT i,j, i_increase,ntries ;
 
+  struct LocGlob lg;
+
+
   /* Assume ok initially: */
   *retval = 0.;
   *epsnegative = 0; 
@@ -275,7 +269,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_Polytype_pt)  (
 
   // Calculate various scalars (Q.B, Q^2, etc)  from the conserved variables:
   
-  Bsq = 
+  lg.Bsq = 
     (*gxx) * (*Bx) * (*Bx) + 
     (*gyy) * (*By) * (*By) + 
     (*gzz) * (*Bz) * (*Bz) + 
@@ -285,25 +279,23 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_Polytype_pt)  (
        (*gyz) * (*By) * (*Bz) );
     
   QdotB = (sx * (*Bx) + sy * (*By) + sz * (*Bz)) ;
-  QdotBsq = QdotB*QdotB ;
-
-  Qtsq = (usx * sx  +  usy * sy  +  usz * sz) ;
+  lg.QdotBsq = QdotB*QdotB ;
 
-  D = dens;
+  lg.Qtsq = (usx * sx  +  usy * sy  +  usz * sz) ;
 
-  half_Bsq = 0.5*Bsq;
+  lg.D = dens;
 
-  t2 = D*D;
-  t4 = QdotBsq*t2;
-  t7 = Bsq*Bsq;
-  t24 = 1/t2;
-  two_Bsq = Bsq + Bsq;
-  t300 = QdotBsq*Bsq*t2;
-  t400 = Qtsq*t2;
+  t2 = lg.D*lg.D;
+  lg.t4 = lg.QdotBsq*t2;
+  lg.t7 = lg.Bsq*lg.Bsq;
+  lg.t24 = 1/t2;
+  lg.two_Bsq = lg.Bsq + lg.Bsq;
+  lg.t300 = lg.QdotBsq*lg.Bsq*t2;
+  lg.t400 = lg.Qtsq*t2;
 
-  s200 = D*gammaeos*(*sc_in);
+  lg.s200 = lg.D*gammaeos*(*sc_in);
   CCTK_REAL g_o_gm1 = (gammaeos/(gammaeos-1.));
-  s100 = g_o_gm1*(*sc_in);
+  lg.s100 = g_o_gm1*(*sc_in);
 
   /* calculate W from last timestep and use for guess */
   vsq = 
@@ -328,7 +320,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_Polytype_pt)  (
 	
   // Always calculate rho from D and gamma so that using D in EOS remains consistent
   //   i.e. you don't get positive values for dP/d(vsq) . 
-  rho0 = D / gamma ;
+  rho0 = lg.D / gamma ;
   CCTK_REAL gm1 = gammaeos-1.;
   rho_gm1 = pow(rho0,gm1);
   p = (*sc_in) * rho_gm1 / gamma; 
@@ -340,20 +332,20 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_Polytype_pt)  (
 
   // Make sure that W is large enough so that v^2 < 1 : 
   i_increase = 0;
-  while( (( W_last*W_last*W_last * ( W_last + 2.*Bsq ) 
-	    - QdotBsq*(2.*W_last + Bsq) ) <= W_last*W_last*(Qtsq-Bsq*Bsq))
+  while( (( W_last*W_last*W_last * ( W_last + 2.*lg.Bsq ) 
+	    - lg.QdotBsq*(2.*W_last + lg.Bsq) ) <= W_last*W_last*(lg.Qtsq-lg.Bsq*lg.Bsq))
 	 && (i_increase < 10) ) {
     W_last *= 10.;
     i_increase++;
   }
 
-  W_for_gnr2 = W_for_gnr2_old = W_last;
-  rho_for_gnr2 = rho_for_gnr2_old = rho0;
+  lg.W_for_gnr2 = lg.W_for_gnr2_old = W_last;
+  lg.rho_for_gnr2 = lg.rho_for_gnr2_old = rho0;
 
   // Calculate W: 
   x_1d[0] = W_last;
 
-  *retval = 1.0*oned_newton_raphson( x_1d, 1, gammaeos, func_W ) ;  
+  *retval = 1.0*oned_newton_raphson( x_1d, 1, gammaeos, &lg, func_W ) ;  
 
   W = x_1d[0];
 	
@@ -369,15 +361,15 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_Polytype_pt)  (
     }
   }
 
-  rho_g = x_rho[0] = rho_for_gnr2; 
+  rho_g = x_rho[0] = lg.rho_for_gnr2; 
 
   ntries = 0;
-  while (  (*retval = gnr2( x_rho, 1, gammaeos, func_rho)) &&  ( ntries++ < 10 )  ) { 
+  while (  (*retval = gnr2( x_rho, 1, gammaeos, &lg, func_rho)) &&  ( ntries++ < 10 )  ) { 
     rho_g *= 10.;
     x_rho[0] = rho_g;
   }
 
-  rho_for_gnr2 = x_rho[0];
+  lg.rho_for_gnr2 = x_rho[0];
 
   if( (*retval != 0) ) {
     *retval = 10;
@@ -385,10 +377,10 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_Polytype_pt)  (
   }
 
   // Calculate v^2 : 
-  rho0 = rho_for_gnr2;
+  rho0 = lg.rho_for_gnr2;
   rho_gm1 = pow(rho0,gm1);
 
-  utsq = (D-rho0)*(D+rho0)/(rho0*rho0);
+  utsq = (lg.D-rho0)*(lg.D+rho0)/(rho0*rho0);
 
   gamma_sq = 1.+utsq;
   gamma = sqrt(gamma_sq);
@@ -421,9 +413,9 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_Polytype_pt)  (
   *w_lorentz = gamma; 
   *pressure = p ; 
 
-  g_o_WBsq = 1./(W+Bsq);
+  g_o_WBsq = 1./(W+lg.Bsq);
   QdB_o_W  = QdotB / W; 
-  *bsq = Bsq * (1.-vsq) + QdB_o_W*QdB_o_W;
+  *bsq = lg.Bsq * (1.-vsq) + QdB_o_W*QdB_o_W;
 
   *velx = g_o_WBsq * ( usx + QdB_o_W*(*Bx) ) ;
   *vely = g_o_WBsq * ( usy + QdB_o_W*(*By) ) ;
@@ -449,47 +441,6 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_Polytype_pt)  (
 }
 
 
-/**********************************************************************/ 
-/****************************************************************************
-   vsq_calc(): 
-    
-      -- evaluate v^2 (spatial, normalized velocity) from 
-            W = \gamma^2 w 
-
-****************************************************************************/
-static CCTK_REAL vsq_calc(CCTK_REAL W)
-{
-	CCTK_REAL Wsq,Xsq,Bsq_W;
-	
-	Wsq = W*W ;
-	Bsq_W = (Bsq + W);
-	Xsq = Bsq_W * Bsq_W;
-
-	return(  ( Wsq * Qtsq  + QdotBsq * (Bsq_W + W)) / (Wsq*Xsq) );
-}
-
-
-/********************************************************************
-
-  x1_of_x0(): 
-           
-    -- calculates v^2 from W  with some physical bounds checking;
-    -- asumes x0 is already physical
-    -- makes v^2 physical  if not;
-
-*********************************************************************/
-
-static CCTK_REAL x1_of_x0(CCTK_REAL x0 ) 
-{
-  CCTK_REAL x1,vsq;
-  CCTK_REAL dv = 1.e-15;
-
-  vsq = fabs(vsq_calc(x0)) ; // guaranteed to be positive 
-
-  return( ( vsq > 1. ) ? (1.0 - dv) : vsq   ); 
-
-}
-
 /********************************************************************
 
   validate_x(): 
@@ -524,10 +475,10 @@ static void validate_x(CCTK_REAL x[2], CCTK_REAL x0[2] )
     -- inspired in part by Num. Rec.'s routine newt();
 
 *****************************************************************/
-static CCTK_INT oned_newton_raphson( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammaeos,
+static CCTK_INT oned_newton_raphson( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammaeos, struct LocGlob *lgp,
 				   void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], 
 						  CCTK_REAL [][1], CCTK_REAL *, 
-						  CCTK_REAL *, CCTK_INT, CCTK_REAL) )
+						  CCTK_REAL *, CCTK_INT, CCTK_REAL, struct LocGlob *) )
 {
   CCTK_REAL f, df, dx[1], x_old[1], resid[1], 
     jac[1][1];
@@ -554,7 +505,7 @@ static CCTK_INT oned_newton_raphson( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammae
   keep_iterating = 1;
   while( keep_iterating ) { 
 
-    (*funcd) (x, dx, resid, jac, &f, &df, n, gammaeos);  /* returns with new dx, f, df */
+    (*funcd) (x, dx, resid, jac, &f, &df, n, gammaeos, lgp);  /* returns with new dx, f, df */
 
     /* Save old values before calculating the new: */
     errx = 0.;
@@ -568,8 +519,8 @@ static CCTK_INT oned_newton_raphson( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammae
 
 //    //METHOD specific:
 //    i_increase = 0;
-//    while( (( x[0]*x[0]*x[0] * ( x[0] + 2.*Bsq ) - 
-//	      QdotBsq*(2.*x[0] + Bsq) ) <= x[0]*x[0]*(Qtsq-Bsq*Bsq))
+//    while( (( x[0]*x[0]*x[0] * ( x[0] + 2.*lgp->Bsq ) - 
+//	      lgp->QdotBsq*(2.*x[0] + lgp->Bsq) ) <= x[0]*x[0]*(lgp->Qtsq-lgp->Bsq*lgp->Bsq))
 //	   && (i_increase < 10) ) {
 //      x[0] -= (1.*i_increase) * dx[0] / 10. ;
 //      i_increase++;
@@ -615,7 +566,7 @@ static CCTK_INT oned_newton_raphson( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammae
   if( (!finite(f)) || (!finite(df)) || (!finite(x[0]))  ) {
 #if(DEBUG_CON2PRIMM)
     fprintf(stderr,"\ngnr not finite, f,df,x_o,x,W_o,W,rho_o,rho = %26.20e %26.20e %26.20e %26.20e %26.20e %26.20e %26.20e %26.20e \n",
-    	    f,df,x[0],x_old[0],W_for_gnr2_old,W_for_gnr2,rho_for_gnr2_old,rho_for_gnr2); fflush(stderr); 
+    	    f,df,x[0],x_old[0],lgp->W_for_gnr2_old,lgp->W_for_gnr2,lgp->rho_for_gnr2_old,lgp->rho_for_gnr2); fflush(stderr); 
 #endif
     return(2);
   }
@@ -652,52 +603,51 @@ static CCTK_INT oned_newton_raphson( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammae
          n    = dimension of x[];
  *********************************************************************************/
 static void func_W(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], 
-		   CCTK_REAL jac[][1], CCTK_REAL *f, CCTK_REAL *df, CCTK_INT n, CCTK_REAL gammaeos)
+		   CCTK_REAL jac[][1], CCTK_REAL *f, CCTK_REAL *df, CCTK_INT n, CCTK_REAL gammaeos, 
+		   struct LocGlob *lgp)
 {
   CCTK_INT retval, ntries;
   CCTK_REAL  W, x_rho[1], rho, rho_g ;
-  CCTK_REAL  t15,   t200,t1000 ;
+  CCTK_REAL  t15,   t200,t1000, drho_dW, rho_gm1 ;
   
   W  = x[0];
-  W_for_gnr2_old = W_for_gnr2;
-  W_for_gnr2 = W;
+  lgp->W_for_gnr2_old = lgp->W_for_gnr2;
+  lgp->W_for_gnr2 = W;
 
   // get rho from NR:
-  rho_g = x_rho[0] = rho_for_gnr2;
+  rho_g = x_rho[0] = lgp->rho_for_gnr2;
   
   ntries = 0;
-  while (  (retval = gnr2( x_rho, 1, gammaeos, func_rho)) &&  ( ntries++ < 10 )  ) { 
+  while (  (retval = gnr2( x_rho, 1, gammaeos, lgp, func_rho)) &&  ( ntries++ < 10 )  ) { 
     rho_g *= 10.;
     x_rho[0] = rho_g;
   }
 
 #if(DEBUG_CON2PRIMM)
   if( x_rho[0] <= 0. ) { 
-    fprintf(stderr,"gnr2 neg rho = %d ,rho_n,rho,rho_o,W,W_o = %26.20e %26.20e %26.20e %26.20e %26.20e \n", retval, x_rho[0], rho_for_gnr2, rho_for_gnr2_old, x[0], W_for_gnr2_old);
+    fprintf(stderr,"gnr2 neg rho = %d ,rho_n,rho,rho_o,W,W_o = %26.20e %26.20e %26.20e %26.20e %26.20e \n", retval, x_rho[0], lgp->rho_for_gnr2, lgp->rho_for_gnr2_old, x[0], lgp->W_for_gnr2_old);
     fflush(stderr);
   }
   
   if( retval ) { 
-    fprintf(stderr,"gnr2 retval = %d ,rho_n,rho,rho_o,W,W_o = %26.20e %26.20e %26.20e %26.20e %26.20e \n", retval, x_rho[0], rho_for_gnr2, rho_for_gnr2_old, x[0], W_for_gnr2_old);
+    fprintf(stderr,"gnr2 retval = %d ,rho_n,rho,rho_o,W,W_o = %26.20e %26.20e %26.20e %26.20e %26.20e \n", retval, x_rho[0], lgp->rho_for_gnr2, lgp->rho_for_gnr2_old, x[0], lgp->W_for_gnr2_old);
     fflush(stderr);
   }
 #endif 
 
-  rho_for_gnr2_old = rho_for_gnr2; 
-  rho = rho_for_gnr2 = x_rho[0];
+  lgp->rho_for_gnr2_old = lgp->rho_for_gnr2; 
+  rho = lgp->rho_for_gnr2 = x_rho[0];
 
   CCTK_REAL gm1 = gammaeos-1.;
 
   rho_gm1 = pow(rho,gm1);
-  drho_dW = -rho*rho/( -rho_gm1*s200 + W*rho);
+  drho_dW = -rho*rho/( -rho_gm1*lgp->s200 + W*rho);
 
-  //  t6 = rho*rho;
-  //  t100 = (D-rho)*(D+rho);  // t2 - t6;
-  t15 = -(D-rho)*(D+rho);  // t6-t2
-  t200 = W + two_Bsq;
+  t15 = -(lgp->D-rho)*(lgp->D+rho);  // t6-t2
+  t200 = W + lgp->two_Bsq;
   t1000 = rho*drho_dW;
-  resid[0] = (t300+(t4+t4+(t400+t15*(t7+(t200)*W))*W)*W)*t24;
-  jac[0][0] = 2*(t4+(t400+t15*t7+(3.0*t15*Bsq+t7*t1000+(t15+t15+t1000*(t200))*W)*W)*W)*t24;
+  resid[0] = (lgp->t300+(lgp->t4+lgp->t4+(lgp->t400+t15*(lgp->t7+(t200)*W))*W)*W)*lgp->t24;
+  jac[0][0] = 2*(lgp->t4+(lgp->t400+t15*lgp->t7+(3.0*t15*lgp->Bsq+lgp->t7*t1000+(t15+t15+t1000*(t200))*W)*W)*W)*lgp->t24;
 
   dx[0] = -resid[0]/jac[0][0];
 
@@ -705,11 +655,11 @@ static void func_W(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[],
   *f = -0.5*(*df);
 
 
-//  fprintf(stdout,"QdotBsq = %28.18e ; \n",QdotBsq );
+//  fprintf(stdout,"QdotBsq = %28.18e ; \n",lgp->QdotBsq );
 //  fprintf(stdout,"Sc      = %28.18e ; \n",Sc     );
-//  fprintf(stdout,"Bsq     = %28.18e ; \n",Bsq     );
-//  fprintf(stdout,"Qtsq    = %28.18e ; \n",Qtsq    );
-//  fprintf(stdout,"Dc      = %28.18e ; \n",D       );
+//  fprintf(stdout,"Bsq     = %28.18e ; \n",lgp->Bsq     );
+//  fprintf(stdout,"Qtsq    = %28.18e ; \n",lgp->Qtsq    );
+//  fprintf(stdout,"Dc      = %28.18e ; \n",lgp->D       );
 //  fprintf(stdout,"drhodW  = %28.18e ; \n",drho_dW  );
 //  fprintf(stdout,"W       = %28.18e ; \n",W       );
 //  fprintf(stdout,"rho     = %28.18e ; \n",rho     );
@@ -730,9 +680,9 @@ static void func_W(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[],
     -- used to calculate rho from W
 
 *****************************************************************/
-static CCTK_INT gnr2( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammaeos, 
+static CCTK_INT gnr2( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammaeos, struct LocGlob *lgp,
 			    void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], 
-					   CCTK_REAL [][1],CCTK_REAL *,CCTK_REAL *,CCTK_INT,CCTK_REAL) )
+					   CCTK_REAL [][1],CCTK_REAL *,CCTK_REAL *,CCTK_INT,CCTK_REAL, struct LocGlob *lgp) )
 {
   CCTK_REAL f, df, dx[1], x_old[1], resid[1], 
     jac[1][1];
@@ -757,7 +707,7 @@ static CCTK_INT gnr2( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammaeos,
   keep_iterating = 1;
   while( keep_iterating ) { 
 
-    (*funcd) (x, dx, resid, jac, &f, &df, n,gammaeos);  /* returns with new dx, f, df */
+    (*funcd) (x, dx, resid, jac, &f, &df, n,gammaeos, lgp);  /* returns with new dx, f, df */
 
     /* Save old values before calculating the new: */
     //-fast  errx = 0.;
@@ -801,7 +751,7 @@ static CCTK_INT gnr2( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammaeos,
   if( (!finite(f)) || (!finite(df)) || (!finite(x[0]))  ) {
 #if(DEBUG_CON2PRIMM)
     fprintf(stderr,"\ngnr2 not finite, f,df,x_o,x,W_o,W,rho_o,rho = %26.20e %26.20e %26.20e %26.20e %26.20e %26.20e %26.20e %26.20e \n",
-    	    f,df,x[0],x_old[0],W_for_gnr2_old,W_for_gnr2,rho_for_gnr2_old,rho_for_gnr2); fflush(stderr); 
+    	    f,df,x[0],x_old[0],lgp->W_for_gnr2_old,lgp->W_for_gnr2,lgp->rho_for_gnr2_old,lgp->rho_for_gnr2); fflush(stderr); 
 #endif
     return(2);
   }
@@ -839,7 +789,8 @@ static CCTK_INT gnr2( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammaeos,
  *********************************************************************************/
 // for the isentropic version:   eq.  (27)
 static void func_rho(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], 
-		     CCTK_REAL jac[][1], CCTK_REAL *f, CCTK_REAL *df, CCTK_INT n, CCTK_REAL gammaeos)
+		     CCTK_REAL jac[][1], CCTK_REAL *f, CCTK_REAL *df, CCTK_INT n, CCTK_REAL gammaeos, 
+		     struct LocGlob *lgp)
 {
 
   CCTK_REAL A, B, C, rho, W, B0;
@@ -848,13 +799,13 @@ static void func_rho(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[],
   CCTK_REAL gm1 = gammaeos-1.;
 
   rho = x[0];
-  W = W_for_gnr2;
+  W = lgp->W_for_gnr2;
 
-  rho_gm1 = t40 = pow(rho,gm1);
+  t40 = pow(rho,gm1);
 
-  resid[0] = (rho*W+(-t40*s100-D)*D);
+  resid[0] = (rho*W+(-t40*lgp->s100-lgp->D)*lgp->D);
   t14 = t40/rho;  // rho^(g-2)
-  jac[0][0] = -t14*s200 + W;
+  jac[0][0] = -t14*lgp->s200 + W;
   //  drho_dW = -rho/jac[0][0];
 
   dx[0] = -resid[0]/jac[0][0];
@@ -862,7 +813,7 @@ static void func_rho(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[],
   *f = -0.5*(*df);
 
   //  fprintf(stdout,"deriv3 %g %g %g %g %g \n",rho,W,resid[0],jac[0][0],dx[0]);
-//  fprintf(stdout,"Dc        := %28.18e ; \n",D);
+//  fprintf(stdout,"Dc        := %28.18e ; \n",lgp->D);
 //  fprintf(stdout,"Sc        := %28.18e ; \n",Sc);
 //  fprintf(stdout,"rho       := %28.18e ; \n",rho);
 //  fprintf(stdout,"W         := %28.18e ; \n",W);
@@ -874,40 +825,6 @@ static void func_rho(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[],
 }
 
 
-/********************************************************************** 
- ********************************************************************** 
-   
- The following routines specify the equation of state.  All routines 
-  above here should be indpendent of EOS.  If the user wishes 
-  to use another equation of state, the below functions must be replaced 
-  by equivalent routines based upon the new EOS. 
-
- **********************************************************************
-**********************************************************************/
-
-/**********************************************************************/
-/********************************************************************** 
-  eos_info():
- 
-      -- returns with all the EOS-related values needed;
- **********************************************************************/
-static CCTK_REAL eos_info(CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL *dpdw, CCTK_REAL *dpdvsq, CCTK_REAL gammaeos)
-{
-  register double ftmp,gtmp;
-
-  ftmp = 1. - vsq;
-  gtmp = sqrt(ftmp);
-
-  CCTK_REAL gam_m1_o_gam = (gammaeos-1.)/gammaeos;
-
-  *dpdw =  gam_m1_o_gam * ftmp ;
-  *dpdvsq =  gam_m1_o_gam * ( 0.5 * D/gtmp  -  W ) ;
-
-  return( gam_m1_o_gam * ( W * ftmp  -  D * gtmp )  );  // p 
-
-}
-
-
 /****************************************************************************** 
              END   
  ******************************************************************************/
diff --git a/src/GRHydro_Con2PrimM_pt.c b/src/GRHydro_Con2PrimM_pt.c
index f5ffd56..cfcf4c0 100644
--- a/src/GRHydro_Con2PrimM_pt.c
+++ b/src/GRHydro_Con2PrimM_pt.c
@@ -71,23 +71,26 @@
 ***************************************************/
 
 /* Local Globals */
-static CCTK_REAL Bsq,QdotBsq,Qtsq,Qdotn,D,half_Bsq ;
+struct LocGlob {
+  CCTK_REAL Bsq,QdotBsq,Qtsq,Qdotn,D,half_Bsq ;
+} ;
 
 // Declarations: 
-static CCTK_REAL vsq_calc(CCTK_REAL W);
+static CCTK_REAL vsq_calc(CCTK_REAL W, struct LocGlob *lgp);
 
-static CCTK_INT twod_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos,
+static CCTK_INT twod_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos, struct LocGlob *lgp,
 				     void (*funcd) (CCTK_REAL [], CCTK_REAL [], 
 						    CCTK_REAL [], CCTK_REAL [][2], 
-						    CCTK_REAL *, CCTK_REAL *, CCTK_REAL) );
+						    CCTK_REAL *, CCTK_REAL *, CCTK_REAL, struct LocGlob *) );
 
 static  void func_vsq( CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], CCTK_REAL [][2], 
-		       CCTK_REAL *f, CCTK_REAL *df, CCTK_REAL gammaeos);
-static CCTK_REAL x1_of_x0(CCTK_REAL x0 ) ;
+		       CCTK_REAL *f, CCTK_REAL *df, CCTK_REAL gammaeos, struct LocGlob *lgp);
+
+static CCTK_REAL x1_of_x0(CCTK_REAL x0, struct LocGlob *lgp ) ;
 
 
 // EOS STUFF:
-static CCTK_REAL eos_info(CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL *dpdw, CCTK_REAL *dpdvsq, CCTK_REAL gammaeos);
+static CCTK_REAL eos_info(CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL *dpdw, CCTK_REAL *dpdvsq, CCTK_REAL gammaeos, struct LocGlob *lgp);
 /* pressure as a function of rho0 and u */
 static CCTK_REAL pressure_rho0_u(CCTK_REAL rho0, CCTK_REAL u, CCTK_REAL gammaeos)
 {
@@ -213,6 +216,8 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt)  (
   CCTK_REAL sqrt_detg = sqrt(detg);
   CCTK_REAL inv_sqrt_detg = 1./sqrt_detg; 
   CCTK_INT i,j, i_increase ;
+  
+  struct LocGlob lg; 
 
   gammaeos = *gamma_eos;
 
@@ -268,7 +273,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt)  (
 
   // Calculate various scalars (Q.B, Q^2, etc)  from the conserved variables:
   
-  Bsq = 
+  lg.Bsq = 
     (*gxx) * (*Bx) * (*Bx) + 
     (*gyy) * (*By) * (*By) + 
     (*gzz) * (*Bz) * (*Bz) + 
@@ -278,15 +283,15 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt)  (
        (*gyz) * (*By) * (*Bz) );
     
   QdotB = (sx * (*Bx) + sy * (*By) + sz * (*Bz)) ;
-  QdotBsq = QdotB*QdotB ;
+  lg.QdotBsq = QdotB*QdotB ;
 
-  Qdotn = -(tau + dens) ;
+  lg.Qdotn = -(tau + dens) ;
 
-  Qtsq = (usx * sx  +  usy * sy  +  usz * sz) ;
+  lg.Qtsq = (usx * sx  +  usy * sy  +  usz * sz) ;
 
-  D = dens;
+  lg.D = dens;
 
-  half_Bsq = 0.5*Bsq;
+  lg.half_Bsq = 0.5*lg.Bsq;
 
   /* calculate W from last timestep and use for guess */
   vsq = 
@@ -311,7 +316,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt)  (
 	
   // Always calculate rho from D and gamma so that using D in EOS remains consistent
   //   i.e. you don't get positive values for dP/d(vsq) . 
-  rho0 = D / gamma ;
+  rho0 = lg.D / gamma ;
   u = (*epsilon) * rho0;
   p = pressure_rho0_u(rho0,u,gammaeos) ;  // EOS
   w = rho0 + u + p ;
@@ -321,8 +326,8 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt)  (
 
   // Make sure that W is large enough so that v^2 < 1 : 
   i_increase = 0;
-  while( (( W_last*W_last*W_last * ( W_last + 2.*Bsq ) 
-	    - QdotBsq*(2.*W_last + Bsq) ) <= W_last*W_last*(Qtsq-Bsq*Bsq))
+  while( (( W_last*W_last*W_last * ( W_last + 2.*lg.Bsq ) 
+	    - lg.QdotBsq*(2.*W_last + lg.Bsq) ) <= W_last*W_last*(lg.Qtsq-lg.Bsq*lg.Bsq))
 	 && (i_increase < 10) ) {
     W_last *= 10.;
     i_increase++;
@@ -330,8 +335,8 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt)  (
   
   // Calculate W and vsq: 
   x_2d[0] =  fabs( W_last );
-  x_2d[1] = x1_of_x0( W_last ) ;
-  *retval = 1.0*twod_newton_raphson( x_2d, gammaeos, func_vsq ) ;  
+  x_2d[1] = x1_of_x0( W_last, &lg ) ;
+  *retval = 1.0*twod_newton_raphson( x_2d, gammaeos, &lg, func_vsq ) ;  
 
   W = x_2d[0];
   vsq = x_2d[1];
@@ -357,7 +362,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt)  (
   // Recover the primitive variables from the scalars and conserved variables:
   gtmp = sqrt(1. - vsq);
   gamma = 1./gtmp ;
-  rho0 = D * gtmp;
+  rho0 = lg.D * gtmp;
 
   w = W * (1. - vsq) ;
   p = pressure_rho0_w(rho0,w,gammaeos) ;  // EOS
@@ -375,9 +380,9 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt)  (
   *w_lorentz = gamma; 
   *pressure = p ; 
 
-  g_o_WBsq = 1./(W+Bsq);
+  g_o_WBsq = 1./(W+lg.Bsq);
   QdB_o_W  = QdotB / W; 
-  *bsq = Bsq * (1.-vsq) + QdB_o_W*QdB_o_W;
+  *bsq = lg.Bsq * (1.-vsq) + QdB_o_W*QdB_o_W;
 
   *velx = g_o_WBsq * ( usx + QdB_o_W*(*Bx) ) ;
   *vely = g_o_WBsq * ( usy + QdB_o_W*(*By) ) ;
@@ -411,15 +416,15 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt)  (
             W = \gamma^2 w 
 
 ****************************************************************************/
-static CCTK_REAL vsq_calc(CCTK_REAL W)
+static CCTK_REAL vsq_calc(CCTK_REAL W, struct LocGlob *lgp)
 {
 	CCTK_REAL Wsq,Xsq,Bsq_W;
 	
 	Wsq = W*W ;
-	Bsq_W = (Bsq + W);
+	Bsq_W = (lgp->Bsq + W);
 	Xsq = Bsq_W * Bsq_W;
 
-	return(  ( Wsq * Qtsq  + QdotBsq * (Bsq_W + W)) / (Wsq*Xsq) );
+	return(  ( Wsq * lgp->Qtsq  + lgp->QdotBsq * (Bsq_W + W)) / (Wsq*Xsq) );
 }
 
 
@@ -433,12 +438,12 @@ static CCTK_REAL vsq_calc(CCTK_REAL W)
 
 *********************************************************************/
 
-static CCTK_REAL x1_of_x0(CCTK_REAL x0 ) 
+static CCTK_REAL x1_of_x0(CCTK_REAL x0, struct LocGlob *lgp ) 
 {
   CCTK_REAL x1,vsq;
   CCTK_REAL dv = 1.e-15;
 
-  vsq = fabs(vsq_calc(x0)) ; // guaranteed to be positive 
+  vsq = fabs(vsq_calc(x0,lgp)) ; // guaranteed to be positive 
 
   return( ( vsq > 1. ) ? (1.0 - dv) : vsq   ); 
 
@@ -478,10 +483,10 @@ static void validate_x(CCTK_REAL x[2], CCTK_REAL x0[2] )
     -- inspired in part by Num. Rec.'s routine newt();
 
 *****************************************************************/
-static CCTK_INT twod_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos, 
+static CCTK_INT twod_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos, struct LocGlob *lgp, 
 			    void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], 
 					   CCTK_REAL [][2], CCTK_REAL *, 
-					   CCTK_REAL *, CCTK_REAL) )
+					   CCTK_REAL *, CCTK_REAL, struct LocGlob *) )
 {
   CCTK_REAL f, df, dx[2], x_old[2];
   CCTK_REAL resid[2], jac[2][2];
@@ -507,7 +512,7 @@ static CCTK_INT twod_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos,
   keep_iterating = 1;
   while( keep_iterating ) { 
 
-    (*funcd) (x, dx, resid, jac, &f, &df, gammaeos);  /* returns with new dx, f, df */
+    (*funcd) (x, dx, resid, jac, &f, &df, gammaeos, lgp);  /* returns with new dx, f, df */
       
 
     /* Save old values before calculating the new: */
@@ -599,7 +604,7 @@ static CCTK_INT twod_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos,
  *********************************************************************************/
 
 static void func_vsq(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], 
-		     CCTK_REAL jac[][2], CCTK_REAL *f, CCTK_REAL *df, CCTK_REAL gammaeos)
+		     CCTK_REAL jac[][2], CCTK_REAL *f, CCTK_REAL *df, CCTK_REAL gammaeos, struct LocGlob *lgp)
 {
 
   
@@ -612,22 +617,22 @@ static void func_vsq(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[],
   
   Wsq = W*W;
   
-  p_tmp = eos_info(W, vsq, &dPdW, &dPdvsq, gammaeos);
+  p_tmp = eos_info(W, vsq, &dPdW, &dPdvsq, gammaeos, lgp);
 
   // These expressions were calculated using Mathematica, but made into efficient 
   // code using Maple.  Since we know the analytic form of the equations, we can 
   // explicitly calculate the Newton-Raphson step: 
 
-  t2 = -half_Bsq+dPdvsq;
-  t3 = Bsq+W;
+  t2 = -lgp->half_Bsq+dPdvsq;
+  t3 = lgp->Bsq+W;
   t4 = t3*t3;
   t23 = 1/W;
-  t16 = QdotBsq*t23*t23;
-  t11 = Qtsq-vsq*t4+t16*(Bsq+W+W);
-  t18 = -Qdotn-half_Bsq*(1.0+vsq)+0.5*t16-W+p_tmp;
+  t16 = lgp->QdotBsq*t23*t23;
+  t11 = lgp->Qtsq-vsq*t4+t16*(lgp->Bsq+W+W);
+  t18 = -lgp->Qdotn-lgp->half_Bsq*(1.0+vsq)+0.5*t16-W+p_tmp;
   t24 = t16*t23;
   t25 = -1.0+dPdW-t24;
-  t35 = t25*t3+(Bsq-2.0*dPdvsq)*(t16+vsq*W)*t23;
+  t35 = t25*t3+(lgp->Bsq-2.0*dPdvsq)*(t16+vsq*W)*t23;
   //  t21 = 1/t3;
   //  t36 = 1/t35;
   t21 = 1/(t3*t35);
@@ -666,7 +671,7 @@ static void func_vsq(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[],
  
       -- returns with all the EOS-related values needed;
  **********************************************************************/
-static CCTK_REAL eos_info(CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL *dpdw, CCTK_REAL *dpdvsq, CCTK_REAL gammaeos)
+static CCTK_REAL eos_info(CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL *dpdw, CCTK_REAL *dpdvsq, CCTK_REAL gammaeos, struct LocGlob *lgp)
 {
   register double ftmp,gtmp;
 
@@ -676,9 +681,9 @@ static CCTK_REAL eos_info(CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL *dpdw, CCTK_REAL
   double gam_m1_o_gam = ((gammaeos-1.)/gammaeos);
 
   *dpdw =  gam_m1_o_gam * ftmp ;
-  *dpdvsq =  gam_m1_o_gam * ( 0.5 * D/gtmp  -  W ) ;
+  *dpdvsq =  gam_m1_o_gam * ( 0.5 * lgp->D/gtmp  -  W ) ;
 
-  return( gam_m1_o_gam * ( W * ftmp  -  D * gtmp )  );  // p 
+  return( gam_m1_o_gam * ( W * ftmp  -  lgp->D * gtmp )  );  // p 
 
 }