diff options
Diffstat (limited to 'src/GRHydro_Con2PrimM_pt.c')
-rw-r--r-- | src/GRHydro_Con2PrimM_pt.c | 100 |
1 files changed, 50 insertions, 50 deletions
diff --git a/src/GRHydro_Con2PrimM_pt.c b/src/GRHydro_Con2PrimM_pt.c index 865defa..b980887 100644 --- a/src/GRHydro_Con2PrimM_pt.c +++ b/src/GRHydro_Con2PrimM_pt.c @@ -61,9 +61,9 @@ #define NEWT_TOL2 (1.0e-15) /* TOL of new 1D^*_{v^2} gnr2 method */ #define MIN_NEWT_TOL2 (1.0e-10) /* TOL of new 1D^*_{v^2} gnr2 method */ -#define W_TOO_BIG (1.e20) /* \gamma^2 (\rho_0 + u + p) is assumed +#define W_TOO_BIG (1.e20) /* \gamma^2 (\rho_0 + u + p) is assumed to always be smaller than this. This - is used to detect solver failures */ + is used to detect solver failures */ #define FAIL_VAL (1.e30) /* Generic value to which we set variables when a problem arises */ @@ -80,12 +80,12 @@ struct LocGlob { static CCTK_REAL vsq_calc(CCTK_REAL W, struct LocGlob *lgp); 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, struct LocGlob *) ); + void (*funcd) (CCTK_REAL [], CCTK_REAL [], + CCTK_REAL [], CCTK_REAL [][2], + 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, struct LocGlob *lgp); + CCTK_REAL *f, CCTK_REAL *df, CCTK_REAL gammaeos, struct LocGlob *lgp); static CCTK_REAL x1_of_x0(CCTK_REAL x0, struct LocGlob *lgp ) ; @@ -177,11 +177,11 @@ RETURN VALUE: of retval = (i*100 + j) where given tolerances; 2 -> Newton-Raphson procedure encountered a numerical divergence (occurrence of "nan" or "+/-inf" ; - + j = 0 -> success 1 -> failure: some sort of failure in Newton-Raphson; 2 -> failure: unphysical vsq = v^2 value at initial guess; - 3 -> failure: W<0 or W>W_TOO_BIG + 3 -> failure: W<0 or W>W_TOO_BIG 4 -> failure: v^2 > 1 ( used to be 5 -> failure: rho,uu <= 0 but now sets epsnegative to non-zero ) @@ -235,37 +235,37 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( #if(DEBUG_CON2PRIMM) fprintf(stdout," *dens = %26.16e \n", *dens_in ); - fprintf(stdout," *sx = %26.16e \n", *sx_in ); - fprintf(stdout," *sy = %26.16e \n", *sy_in ); - fprintf(stdout," *sz = %26.16e \n", *sz_in ); - fprintf(stdout," *tau = %26.16e \n", *tau_in ); - fprintf(stdout," *Bconsx = %26.16e \n", *Bconsx_in ); - fprintf(stdout," *Bconsy = %26.16e \n", *Bconsy_in ); + fprintf(stdout," *sx = %26.16e \n", *sx_in ); + fprintf(stdout," *sy = %26.16e \n", *sy_in ); + fprintf(stdout," *sz = %26.16e \n", *sz_in ); + fprintf(stdout," *tau = %26.16e \n", *tau_in ); + fprintf(stdout," *Bconsx = %26.16e \n", *Bconsx_in ); + fprintf(stdout," *Bconsy = %26.16e \n", *Bconsy_in ); fprintf(stdout," *Bconsz = %26.16e \n", *Bconsz_in ); - fprintf(stdout," *rho = %26.16e \n", *rho ); - fprintf(stdout," *velx = %26.16e \n", *velx ); - fprintf(stdout," *vely = %26.16e \n", *vely ); - fprintf(stdout," *velz = %26.16e \n", *velz ); + fprintf(stdout," *rho = %26.16e \n", *rho ); + fprintf(stdout," *velx = %26.16e \n", *velx ); + fprintf(stdout," *vely = %26.16e \n", *vely ); + fprintf(stdout," *velz = %26.16e \n", *velz ); fprintf(stdout," *epsilon = %26.16e \n", *epsilon ); fprintf(stdout," *pressure = %26.16e \n", *pressure ); - fprintf(stdout," *Bx = %26.16e \n", *Bx ); - fprintf(stdout," *By = %26.16e \n", *By ); - fprintf(stdout," *Bz = %26.16e \n", *Bz ); - fprintf(stdout," *bsq = %26.16e \n", *bsq ); + fprintf(stdout," *Bx = %26.16e \n", *Bx ); + fprintf(stdout," *By = %26.16e \n", *By ); + fprintf(stdout," *Bz = %26.16e \n", *Bz ); + fprintf(stdout," *bsq = %26.16e \n", *bsq ); fprintf(stdout," *w_lorentz = %26.16e \n", *w_lorentz ); - fprintf(stdout," *gxx = %26.16e \n", *gxx ); - fprintf(stdout," *gxy = %26.16e \n", *gxy ); - fprintf(stdout," *gxz = %26.16e \n", *gxz ); - fprintf(stdout," *gyy = %26.16e \n", *gyy ); - fprintf(stdout," *gyz = %26.16e \n", *gyz ); - fprintf(stdout," *gzz = %26.16e \n", *gzz ); - fprintf(stdout," *uxx = %26.16e \n", *uxx ); - fprintf(stdout," *uxy = %26.16e \n", *uxy ); - fprintf(stdout," *uxz = %26.16e \n", *uxz ); - fprintf(stdout," *uyy = %26.16e \n", *uyy ); - fprintf(stdout," *uyz = %26.16e \n", *uyz ); - fprintf(stdout," *uzz = %26.16e \n", *uzz ); - fprintf(stdout," *det = %26.16e \n", *det ); + fprintf(stdout," *gxx = %26.16e \n", *gxx ); + fprintf(stdout," *gxy = %26.16e \n", *gxy ); + fprintf(stdout," *gxz = %26.16e \n", *gxz ); + fprintf(stdout," *gyy = %26.16e \n", *gyy ); + fprintf(stdout," *gyz = %26.16e \n", *gyz ); + fprintf(stdout," *gzz = %26.16e \n", *gzz ); + fprintf(stdout," *uxx = %26.16e \n", *uxx ); + fprintf(stdout," *uxy = %26.16e \n", *uxy ); + fprintf(stdout," *uxz = %26.16e \n", *uxz ); + fprintf(stdout," *uyy = %26.16e \n", *uyy ); + fprintf(stdout," *uyz = %26.16e \n", *uyz ); + fprintf(stdout," *uzz = %26.16e \n", *uzz ); + fprintf(stdout," *det = %26.16e \n", *det ); fprintf(stdout," *epsnegative = %10d \n", *epsnegative ); fprintf(stdout," *retval = %26.16e \n", *retval ); fflush(stdout); @@ -328,7 +328,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( gammasq = 1. / (1. - vsq); gamma = sqrt(gammasq); - + // 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 = lg.D / gamma ; @@ -342,8 +342,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.*lg.Bsq ) - - lg.QdotBsq*(2.*W_last + lg.Bsq) ) <= W_last*W_last*(lg.Qtsq-lg.Bsq*lg.Bsq)) - && (i_increase < 10) ) { + - 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++; } @@ -355,7 +355,7 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( W = x_2d[0]; vsq = x_2d[1]; - + /* Problem with solver, so return denoting error before doing anything further */ if( ((*retval) != 0.) || (W == FAIL_VAL) ) { *retval = *retval*100.+1.; @@ -441,13 +441,13 @@ void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_pt) ( ****************************************************************************/ static CCTK_REAL vsq_calc(CCTK_REAL W, struct LocGlob *lgp) { - CCTK_REAL Wsq,Xsq,Bsq_W; - - Wsq = W*W ; - Bsq_W = (lgp->Bsq + W); - Xsq = Bsq_W * Bsq_W; + CCTK_REAL Wsq,Xsq,Bsq_W; + + Wsq = W*W ; + Bsq_W = (lgp->Bsq + W); + Xsq = Bsq_W * Bsq_W; - return( ( Wsq * lgp->Qtsq + lgp->QdotBsq * (Bsq_W + W)) / (Wsq*Xsq) ); + return( ( Wsq * lgp->Qtsq + lgp->QdotBsq * (Bsq_W + W)) / (Wsq*Xsq) ); } @@ -509,9 +509,9 @@ static void validate_x(CCTK_REAL x[2], CCTK_REAL x0[2] ) *****************************************************************/ 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, struct LocGlob *) ) + void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], + CCTK_REAL [][2], CCTK_REAL *, + CCTK_REAL *, CCTK_REAL, struct LocGlob *) ) { CCTK_REAL f, df, dx[2], x_old[2]; CCTK_REAL resid[2], jac[2][2]; @@ -580,7 +580,7 @@ static CCTK_INT twod_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos, struct L if( doing_extra == 1 ) i_extra++ ; if( ((fabs(errx) <= NEWT_TOL)&&(doing_extra == 0)) - || (i_extra > EXTRA_NEWT_ITER) || (n_iter >= (MAX_NEWT_ITER-1)) ) { + || (i_extra > EXTRA_NEWT_ITER) || (n_iter >= (MAX_NEWT_ITER-1)) ) { keep_iterating = 0; } @@ -629,7 +629,7 @@ static CCTK_INT twod_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos, struct L *********************************************************************************/ 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, struct LocGlob *lgp) + CCTK_REAL jac[][2], CCTK_REAL *f, CCTK_REAL *df, CCTK_REAL gammaeos, struct LocGlob *lgp) { |