aboutsummaryrefslogtreecommitdiff
path: root/src/GRHydro_Con2PrimM_pt.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/GRHydro_Con2PrimM_pt.c')
-rw-r--r--src/GRHydro_Con2PrimM_pt.c85
1 files changed, 45 insertions, 40 deletions
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
}
generated by cgit v1.2.3 (git 2.39.1) at 2024-07-25 16:55:10 +0000