diff options
author | bmundim <bmundim@c83d129a-5a75-4d5a-9c4d-ed3a5855bf45> | 2011-01-06 01:13:02 +0000 |
---|---|---|
committer | bmundim <bmundim@c83d129a-5a75-4d5a-9c4d-ed3a5855bf45> | 2011-01-06 01:13:02 +0000 |
commit | f2cdc7f3a147f13dd3f4861c101f995f8cbfbd05 (patch) | |
tree | b5ea1d82267b2454d6afd70afb4cd3ee962c8d8d /src/GRHydro_Con2PrimM_polytype_pt.c | |
parent | e684ca4eded7391a7a03d5361ed2088b47f38c4f (diff) |
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
Diffstat (limited to 'src/GRHydro_Con2PrimM_polytype_pt.c')
-rw-r--r-- | src/GRHydro_Con2PrimM_polytype_pt.c | 249 |
1 files changed, 83 insertions, 166 deletions
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 ******************************************************************************/ |