[Commits] [svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 208)
bcmsma at astro.rit.edu
bcmsma at astro.rit.edu
Wed Jan 5 19:13:02 CST 2011
User: bmundim
Date: 2011/01/05 07:13 PM
Modified:
/trunk/src/
GRHydro_Con2PrimM.F90, GRHydro_Con2PrimM_polytype_pt.c, GRHydro_Con2PrimM_pt.c
Log:
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.
File Changes:
Directory: /trunk/src/
======================
File [modified]: GRHydro_Con2PrimM.F90
Delta lines: +7 -5
===================================================================
--- trunk/src/GRHydro_Con2PrimM.F90 2011-01-04 14:40:06 UTC (rev 207)
+++ trunk/src/GRHydro_Con2PrimM.F90 2011-01-06 01:13:02 UTC (rev 208)
@@ -94,7 +94,7 @@
!$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 @@
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 @@
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)
File [modified]: GRHydro_Con2PrimM_polytype_pt.c
Delta lines: +83 -166
===================================================================
--- trunk/src/GRHydro_Con2PrimM_polytype_pt.c 2011-01-04 14:40:06 UTC (rev 207)
+++ trunk/src/GRHydro_Con2PrimM_polytype_pt.c 2011-01-06 01:13:02 UTC (rev 208)
@@ -71,43 +71,33 @@
***************************************************/
/* 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) );
+ 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);
+ CCTK_REAL jac[][1], CCTK_REAL *f, CCTK_REAL *df, CCTK_INT n, CCTK_REAL gammaeos,
+ struct LocGlob *lgp);
-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) ;
-}
-
-
void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_Polytype_pt) (
CCTK_INT *handle, CCTK_REAL *gamma_eos,
CCTK_REAL *dens_in,
@@ -220,8 +210,12 @@
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 @@
// 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 @@
(*gyz) * (*By) * (*Bz) );
QdotB = (sx * (*Bx) + sy * (*By) + sz * (*Bz)) ;
- QdotBsq = QdotB*QdotB ;
+ lg.QdotBsq = QdotB*QdotB ;
- 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;
+ 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;
- t2 = D*D;
- t4 = QdotBsq*t2;
- t7 = Bsq*Bsq;
- t24 = 1/t2;
- two_Bsq = Bsq + Bsq;
- t300 = QdotBsq*Bsq*t2;
- t400 = 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 @@
// 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 @@
// 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 @@
}
}
- 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 @@
}
// 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 @@
*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,49 +441,8 @@
}
-/**********************************************************************/
-/****************************************************************************
- 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():
-- makes sure that x[0,1] have physical values, based upon
@@ -524,10 +475,10 @@
-- 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 @@
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 @@
// //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 @@
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 @@
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 @@
*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 @@
-- 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 @@
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 @@
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 @@
*********************************************************************************/
// 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 @@
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 @@
*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 @@
}
-/**********************************************************************
- **********************************************************************
-
- 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
******************************************************************************/
File [modified]: GRHydro_Con2PrimM_pt.c
Delta lines: +45 -40
===================================================================
--- trunk/src/GRHydro_Con2PrimM_pt.c 2011-01-04 14:40:06 UTC (rev 207)
+++ trunk/src/GRHydro_Con2PrimM_pt.c 2011-01-06 01:13:02 UTC (rev 208)
@@ -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 @@
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 @@
// 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 @@
(*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 @@
// 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 @@
// 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 @@
// 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 @@
// 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 @@
*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 @@
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 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 @@
-- 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 @@
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 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 @@
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 @@
-- 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 @@
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
}
More information about the Commits
mailing list