diff --git a/Makefile b/Makefile index d5b54a6..224f48e 100644 --- a/Makefile +++ b/Makefile @@ -20,7 +20,7 @@ endif LIBS=-lm -lgsl -lgslcblas -lsiloh5 -lfftw3 -lrt -lhdf5_serial RM=/bin/rm -OBJS = mpi.o u2prad.o magn.o silo.o postproc.o fileop.o misc.o physics.o finite.o problem.o metric.o relele.o rad.o opacities.o u2p.o frames.o p2u.o nonthermal.o +OBJS = mpi.o u2prad.o magn.o silo.o postproc.o fileop.o misc.o physics.o finite.o problem.o metric.o relele.o rad.o opacities.o u2p.o u2p_ff.o frames.o p2u.o nonthermal.o all: ko ana avg outavg phisli thsli phiavg regrid dumps2hdf5 diff --git a/PROBLEMS/MONOPOLE_2D/bc.c b/PROBLEMS/MONOPOLE_2D/bc.c index 01e15c6..ec70a32 100755 --- a/PROBLEMS/MONOPOLE_2D/bc.c +++ b/PROBLEMS/MONOPOLE_2D/bc.c @@ -112,8 +112,7 @@ pick_G(ix,iy,iz,GG); //in 3D polar cells overwritten with #define CORRECT_POLARAXIS_3D if(BCtype==YBCLO) //upper spin axis { - - iiy=-iy-1; + iiy=-1*iy-1; iiz=iz; iix=ix; gdet_src=get_g(g,3,4,iix,iiy,iiz); @@ -121,7 +120,11 @@ if(BCtype==YBCLO) //upper spin axis for(iv=0;ivgg,geom2->GG); - pp2[2]=ucon[1]; - pp2[3]=ucon[2]; - pp2[4]=ucon[3]; - + pp2[VX]=ucon[1]; + pp2[VY]=ucon[2]; + pp2[VZ]=ucon[3]; + #ifdef MAGNFIELD //coming back to primitive B^i @@ -93,7 +93,26 @@ trans_pmhd_coco(ldouble *ppin, ldouble *ppout, int CO1,int CO2, ldouble *xxvec, { pp2[B1+i] = Bcon[1+i]; } + +#ifdef FORCEFREE + pp2[UUFF]=pp1[UUFF]; + + ldouble uconff[4]; + uconff[0]=0.; + uconff[1]=pp1[VXFF]; + uconff[2]=pp1[VYFF]; + uconff[3]=pp1[VZFF]; + conv_vels(uconff,uconff,VELPRIM,VEL4,geom1->gg,geom1->GG); + trans2_coco(xxvec,uconff,uconff,CO1,CO2); + conv_vels(uconff,uconff,VEL4,VELPRIM,geom2->gg,geom2->GG); + pp2[VXFF]=uconff[1]; + pp2[VYFF]=uconff[2]; + pp2[VZFF]=uconff[3]; +#endif + #endif + + } for(i=0;igg,geom1->GG); + trans2_coco_precompute(uconff,uconff,geom1->ix,geom1->iy,geom1->iz,which); + conv_vels(uconff,uconff,VEL4,VELPRIM,geom2->gg,geom2->GG); + pp2[VXFF]=uconff[1]; + pp2[VYFF]=uconff[2]; + pp2[VZFF]=uconff[3]; +#endif + #endif } diff --git a/ko.c b/ko.c index 1c3ef7d..6dab9dc 100644 --- a/ko.c +++ b/ko.c @@ -114,10 +114,11 @@ main(int argc, char **argv) //print scalings GU->CGS #ifdef FORCEFREE if(PROCID==0) printf("NVHD %d NVMHD %d NV %d\n",NVHD,NVMHD,NV); - if(PROCID==0) printf("%d %d %d %d %d %d %d %d %d %d %d %d\n",RHO,UU,VX,VY,VZ,ENTR,B1,B2,B3,VXFF,VYFF,VZFF); + if(PROCID==0) printf("%d %d %d %d %d %d %d %d %d %d %d %d %d\n", + RHO,UU,VX,VY,VZ,ENTR,B1,B2,B3,UUFF,VXFF,VYFF,VZFF); #endif + if(PROCID==0) print_scalings(); - //exit(1); //************** //tests @@ -245,6 +246,9 @@ main(int argc, char **argv) } calc_BfromA(p,1); +#ifdef FORCEFREE + update_ffprims(); // make force-free primitives consistent +#endif //exchange magn. field calculated in domain mpi_exchangedata(); calc_avgs_throughout(); diff --git a/ko.h b/ko.h index b05c011..f910615 100644 --- a/ko.h +++ b/ko.h @@ -897,23 +897,6 @@ int trans_prad_coco(ldouble *ppin, ldouble *ppout, int CO1,int CO2, ldouble *xxv //deprecated int prad_ff2lab(ldouble *pp1, ldouble *pp2, void* ggg); -//int prad_lab2ff(ldouble *pp1, ldouble *pp2, void *ggg); -//int prad_on2lab(ldouble *pp1, ldouble *pp2, void* ggg); -//int prad_lab2on(ldouble *pp1, ldouble *pp2, void *ggg); -//int prad_ff2zamo(ldouble *pp1, ldouble *pp2, ldouble gg[][5], ldouble GG[][5], ldouble eup[][4]); -//int prad_zamo2ff(ldouble *pp1, ldouble *pp2, ldouble gg[][5], ldouble GG[][5], ldouble eup[][4]); -//int boost2_zamo2ff(ldouble* A1,ldouble* A2,ldouble *pp,ldouble gg[][5],ldouble GG[][5],ldouble eup[][4]); -//int boost2_ff2zamo(ldouble A1[4],ldouble A2[4],ldouble *pp,ldouble gg[][5],ldouble GG[][5],ldouble eup[][4]); -//int boost22_zamo2ff(ldouble T1[][4],ldouble T2[][4],ldouble *pp,ldouble gg[][5],ldouble GG[][5],ldouble eup[][4]); -//int boost22_ff2zamo(ldouble T1[][4],ldouble T2[][4],ldouble *pp,ldouble gg[][5],ldouble GG[][5],ldouble eup[][4]); -//int trans22_zamo2lab(ldouble T1[][4],ldouble T2[][4],ldouble elo[][4]); -//int trans22_lab2zamo(ldouble T1[][4],ldouble T2[][4],ldouble eup[][4]); -//int trans2_lab2zamo(ldouble *u1,ldouble *u2,ldouble eup[][4]); -//int trans2_zamo2lab(ldouble *u1,ldouble *u2,ldouble elo[][4]); -//int trans22_on2cc(ldouble T1[][4],ldouble T2[][4],ldouble tlo[][4]); -//int trans22_cc2on(ldouble T1[][4],ldouble T2[][4],ldouble tup[][4]); -//int trans2_cc2on(ldouble *u1,ldouble *u2,ldouble tup[][4]); -//int trans2_on2cc(ldouble *u1,ldouble *u2,ldouble tlo[][4]); int calc_Lorentz_lab2ff(ldouble *pp,ldouble gg[][5],ldouble GG[][5],ldouble L[][4]); int calc_Lorentz_lab2ff_4vel(ldouble *pp, ldouble gg[][5], ldouble GG[][5], ldouble L[][4], ldouble ucon[4], ldouble ucov[4]); @@ -1024,8 +1007,10 @@ int u2p_solver(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose); int u2p_solver_nonrel(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose); int u2p_solver_Wp(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose); int u2p_solver_W(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose); -int u2p_solver_ff(ldouble *uu, ldouble *pp, void *ggg,int verbose); int u2p_solver_Bonly(ldouble *uu, ldouble *pp, void *ggg); +int u2p_solver_ff(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose); +int update_ffprims(); +int update_ffprims_cell(ldouble *pp, ldouble *uu, void *gg); int count_entropy(int *n, int *n2); int copy_entropycount(); diff --git a/misc.c b/misc.c index 79f5c96..f78b6fc 100644 --- a/misc.c +++ b/misc.c @@ -1476,6 +1476,7 @@ print_primitives(ldouble *p) printf("B^2 = %.15e\n",p[B2]); printf("B^3 = %.15e\n",p[B3]); #ifdef FORCEFREE + printf("uu (ff) = %.15e\n", p[UUFF]); printf("u^1 (ff) = %.15e\n",p[VXFF]); printf("u^2 (ff) = %.15e\n",p[VYFF]); printf("u^3 (ff) = %.15e\n",p[VZFF]); diff --git a/mnemonics.h b/mnemonics.h index e221ff7..90af3df 100644 --- a/mnemonics.h +++ b/mnemonics.h @@ -35,9 +35,10 @@ #define B3 (NVHD+2) #ifdef FORCEFREE -#define VXFF (NVHD+3) -#define VYFF (NVHD+4) -#define VZFF (NVHD+5) +#define UUFF (NVHD+3) +#define VXFF (NVHD+4) +#define VYFF (NVHD+5) +#define VZFF (NVHD+6) #endif #define EE (NVMHD) diff --git a/p2u.c b/p2u.c index 5fb321a..3fd7bb5 100644 --- a/p2u.c +++ b/p2u.c @@ -96,24 +96,13 @@ p2u_mhd(ldouble *pp, ldouble *uu, void *ggg) ldouble Ttr =eta*ucon[0]*ucov[1] - bcon[0]*bcov[1]; ldouble Ttth =eta*ucon[0]*ucov[2] - bcon[0]*bcov[2]; ldouble Ttph =eta*ucon[0]*ucov[3] - bcon[0]*bcov[3]; - /* -#ifdef FORCEFREE - // ANDREW TODO -- we should get rid of this part - - //ldouble Tttt_ff = bsq*ucon[0]*ucov[0] - bcon[0]*bcov[0] + 0.5*bsq + rho*ucon[0]; - Ttr = bsq*ucon[0]*ucov[1] - bcon[0]*bcov[1]; - Ttth = bsq*ucon[0]*ucov[2] - bcon[0]*bcov[2]; - Ttph = bsq*ucon[0]*ucov[3] - bcon[0]*bcov[3]; -#endif - */ uu[RHO]=gdetu*rhout; uu[UU]=gdetu*Tttt; uu[VX]=gdetu*Ttr; uu[VY]=gdetu*Ttth; uu[VZ]=gdetu*Ttph; uu[ENTR]=gdetu*Sut; - #ifdef EVOLVEELECTRONS ldouble Seut=pp[ENTRE]*ut; @@ -125,9 +114,9 @@ p2u_mhd(ldouble *pp, ldouble *uu, void *ggg) #ifdef RELELECTRONS int ib; for(ib=0;ib0) // Save to avg arrays #ifdef DTAVG //Dont save every step @@ -310,7 +310,7 @@ solve_the_problem(ldouble tstart, char* folder) calc_u2p(0,1); // Count number of entropy inversions: ENTROPYFLAG, ENTROPYFLAG2 - count_entropy(&nentr[4],&nentr2[4]); + //count_entropy(&nentr[4],&nentr2[4]); // Special treatment near axis (or inner surface) do_correct(); @@ -333,7 +333,7 @@ solve_the_problem(ldouble tstart, char* folder) global_expdt=dt; // Count number of entropy inversions: ENTROPYFLAG, ENTROPYFLAG2 - count_entropy(&nentr[5],&nentr2[5]); + //count_entropy(&nentr[5],&nentr2[5]); // Calculate 2nd explicit deriv // F(U(2)) in *dut2; @@ -381,7 +381,7 @@ solve_the_problem(ldouble tstart, char* folder) #endif // Count number of entropy inversions: ENTROPYFLAG, ENTROPYFLAG2 - count_entropy(&nentr[6],&nentr2[6]); + //count_entropy(&nentr[6],&nentr2[6]); // ANDREW should the next 3 steps be moved to start of loop? // Special treatment near axis (or inner surface) @@ -397,6 +397,9 @@ solve_the_problem(ldouble tstart, char* folder) // Compute new entropy from rho and uint over the domain update_entropy(); + #ifdef FORCEFREE + update_ffprims(); + #endif } // else if(TIMESTEPPING==RK2IMEX) @@ -416,7 +419,7 @@ solve_the_problem(ldouble tstart, char* folder) copy_u(1.,p,ptm1); // Count number of entropy inversions: ENTROPYFLAG, ENTROPYFLAG2 - count_entropy(&nentr[0],&nentr2[0]); + //count_entropy(&nentr[0],&nentr2[0]); // Special treatment near axis (or inner surface) do_correct(); @@ -435,7 +438,7 @@ solve_the_problem(ldouble tstart, char* folder) op_explicit (t, dt); // Count number of entropy inversions: ENTROPYFLAG, ENTROPYFLAG2 - count_entropy(&nentr[1],&nentr2[1]); + //count_entropy(&nentr[1],&nentr2[1]); copy_entropycount(); // Artifical dynamo (ifdef MIMICDYNAMO) @@ -473,7 +476,7 @@ solve_the_problem(ldouble tstart, char* folder) do_correct(); // Count number of entropy inversions: ENTROPYFLAG, ENTROPYFLAG2 - count_entropy(&nentr[2],&nentr2[2]); + //count_entropy(&nentr[2],&nentr2[2]); // Exchange mpi data for boundaries mpi_exchangedata(); @@ -486,7 +489,7 @@ solve_the_problem(ldouble tstart, char* folder) op_explicit (t,dt); // Count number of entropy inversions: ENTROPYFLAG, ENTROPYFLAG2 - count_entropy(&nentr[3],&nentr2[3]); + //count_entropy(&nentr[3],&nentr2[3]); // Artifical dynamo (ifdef MIMICDYNAMO) apply_dynamo(t,dt); @@ -535,6 +538,10 @@ solve_the_problem(ldouble tstart, char* folder) // Update entropy from rho and uint over the domain update_entropy(); + #ifdef FORCEFREE + update_ffprims(); + #endif + t+=dt; } // else if(TIMESTEPPING==RK2HEUN) @@ -554,7 +561,7 @@ solve_the_problem(ldouble tstart, char* folder) copy_u(1.,p,ptm1); // Count number of entropy inversions: ENTROPYFLAG, ENTROPYFLAG2 - count_entropy(&nentr[0],&nentr2[0]); + //count_entropy(&nentr[0],&nentr2[0]); // Special treatment near axis (or inner surface) do_correct(); @@ -573,7 +580,7 @@ solve_the_problem(ldouble tstart, char* folder) op_explicit (t,.5*dt); // Count number of entropy inversions: ENTROPYFLAG, ENTROPYFLAG2 - count_entropy(&nentr[1],&nentr2[1]); + //count_entropy(&nentr[1],&nentr2[1]); copy_entropycount(); // Artifical dynamo (ifdef MIMICDYNAMO) @@ -611,7 +618,7 @@ solve_the_problem(ldouble tstart, char* folder) do_correct(); // Count number of entropy inversions: ENTROPYFLAG, ENTROPYFLAG2 - count_entropy(&nentr[2],&nentr2[2]); + //count_entropy(&nentr[2],&nentr2[2]); // Exchange MPI data for boundaries mpi_exchangedata(); @@ -624,7 +631,7 @@ solve_the_problem(ldouble tstart, char* folder) op_explicit (t,dt); // Count number of entropy inversions: ENTROPYFLAG, ENTROPYFLAG2 - count_entropy(&nentr[3],&nentr2[3]); + //count_entropy(&nentr[3],&nentr2[3]); // Artifical dynamo (ifdef MIMICDYNAMO) apply_dynamo(t,dt); @@ -673,7 +680,9 @@ solve_the_problem(ldouble tstart, char* folder) // Compute entropy from rho and uint over the domain update_entropy(); - + #ifdef FORCEFREE + update_ffprims(); + #endif t+=dt; } // else if(TIMESTEPPING==RK2) @@ -688,7 +697,7 @@ solve_the_problem(ldouble tstart, char* folder) copyi_u(1.,u,ut0); // Count number of entropy inversions: ENTROPYFLAG, ENTROPYFLAG2 - count_entropy(&nentr[0],&nentr2[0]); + //count_entropy(&nentr[0],&nentr2[0]); // Special treatment near axis (or inner surface) do_correct(); @@ -707,7 +716,8 @@ solve_the_problem(ldouble tstart, char* folder) op_explicit (t,1.*dt); // Count number of entropy inversions: ENTROPYFLAG, ENTROPYFLAG2 - count_entropy(&nentr[1],&nentr2[1]); copy_entropycount(); + //count_entropy(&nentr[1],&nentr2[1]); + copy_entropycount(); // Artifical dynamo (ifdef MIMICDYNAMO) apply_dynamo(t,dt); @@ -751,6 +761,12 @@ solve_the_problem(ldouble tstart, char* folder) heat_electronions_with_state(dt); #endif + // Compute entropy from rho and uint over the domain + update_entropy(); + #ifdef FORCEFREE + update_ffprims(); + #endif + // Update to new time: t+dt t+=dt; } // else if(TIMESTEPPING==RK1) diff --git a/problem.h b/problem.h index ebc2714..7654bea 100644 --- a/problem.h +++ b/problem.h @@ -145,14 +145,26 @@ //141 KEPINF -- INFDISK modified for injecting keplerian //142 PARTIALTDE -- INFDISK modified for partial TDE binding energy distribution non-constant //143 FFTESTS -- Force Free tests +//144 FFTORUSTEST -- test of force-free torus //ANDREW -- I've gone through problems 100-133 and undefined PR_KAPPA where appropriate //If you want to use default calc_opacities_from_state, make sure PR_KAPPA is undefined! //if you are using a problem older than 100 with a different kappa defined in kappa.c //make sure your kappa.c ends with (kappa=(.....)) NOT (return kappa)!! -#define PROBLEM 132 //140 +#define PROBLEM 132 +//#define PROBLEM 144 +#if(PROBLEM==144) + +#define PR_DEFINE "PROBLEMS/FFTORUSTEST/define.h" +#define PR_BC "PROBLEMS/FFTORUSTEST/bc.c" +#define PR_INIT "PROBLEMS/FFTORUSTEST/init.c" +#define PR_KAPPAES "PROBLEMS/FFTORUSTEST/kappaes.c" +#define PR_TOOLS "PROBLEMS/FFTORUSTEST/tools.c" +#define PR_POSTINIT "PROBLEMS/FFTORUSTEST/postinit.c" + +#endif #if(PROBLEM==143) #define PR_DEFINE "PROBLEMS/FFTESTS/define.h" diff --git a/silo.c b/silo.c index 6760b7f..73be334 100644 --- a/silo.c +++ b/silo.c @@ -68,6 +68,7 @@ int fprint_silofile(ldouble time, int num, char* folder, char* prefix) ldouble *rho = (ldouble*)malloc(nx*ny*nz*sizeof(double)); ldouble *entr = (ldouble*)malloc(nx*ny*nz*sizeof(double)); ldouble *uint = (ldouble*)malloc(nx*ny*nz*sizeof(double)); + ldouble *uintff = (ldouble*)malloc(nx*ny*nz*sizeof(double)); ldouble *temp = (ldouble*)malloc(nx*ny*nz*sizeof(double)); ldouble *trad = (ldouble*)malloc(nx*ny*nz*sizeof(double)); ldouble *tradlte = (ldouble*)malloc(nx*ny*nz*sizeof(double)); @@ -83,8 +84,15 @@ int fprint_silofile(ldouble time, int num, char* folder, char* prefix) ldouble *u1 = (ldouble*)malloc(nx*ny*nz*sizeof(double)); ldouble *u2 = (ldouble*)malloc(nx*ny*nz*sizeof(double)); ldouble *u3 = (ldouble*)malloc(nx*ny*nz*sizeof(double)); - - ldouble *lorentz = (ldouble*)malloc(nx*ny*nz*sizeof(double)); + ldouble *lorentz = (ldouble*)malloc(nx*ny*nz*sizeof(double)); + + ldouble *u0_perp = (ldouble*)malloc(nx*ny*nz*sizeof(double)); + ldouble *u1_perp = (ldouble*)malloc(nx*ny*nz*sizeof(double)); + ldouble *u2_perp = (ldouble*)malloc(nx*ny*nz*sizeof(double)); + ldouble *u3_perp = (ldouble*)malloc(nx*ny*nz*sizeof(double)); + ldouble *lorentz_perp = (ldouble*)malloc(nx*ny*nz*sizeof(double)); + ldouble *lorentz_par = (ldouble*)malloc(nx*ny*nz*sizeof(double)); + ldouble *vx = (ldouble*)malloc(nx*ny*nz*sizeof(double)); ldouble *vy = (ldouble*)malloc(nx*ny*nz*sizeof(double)); ldouble *vz = (ldouble*)malloc(nx*ny*nz*sizeof(double)); @@ -740,6 +748,28 @@ int fprint_silofile(ldouble time, int num, char* folder, char* prefix) u2[zonalindex]=vel[2]; u3[zonalindex]=vel[3]; +#ifdef FORCEFREE // must have VELPRIM=VELR + ldouble velff[4]; + velff[0] = 0.; + velff[1]=pp[VXFF]; + velff[2]=pp[VYFF]; + velff[3]=pp[VZFF]; + conv_vels(velff,velff,VELR,VEL4,geomout.gg,geomout.GG); + lorentz_perp[zonalindex]=fabs(velff[0])/sqrt(fabs(geomout.GG[0][0])); + lorentz_par[zonalindex]=1./sqrt(1+pow(lorentz[zonalindex],-2)-pow(lorentz_perp[zonalindex],-2)); + + u0_perp[zonalindex]=velff[0]; + u1_perp[zonalindex]=velff[1]; + u2_perp[zonalindex]=velff[2]; + u3_perp[zonalindex]=velff[3]; + + uintff[zonalindex]=pp[UUFF]; +#else + //TODO actually get || and perp velocities here + lorentz_perp[zonalindex]=-1; + lorentz_par[zonalindex]=-1; +#endif + #ifdef EVOLVEELECTRONS tempi[zonalindex]=tempiloc; //ion temperature tempe[zonalindex]=tempeloc; //electron temperature @@ -1531,6 +1561,32 @@ int fprint_silofile(ldouble time, int num, char* folder, char* prefix) dimensions, ndim, NULL, 0, DB_DOUBLE, DB_ZONECENT, optList); +#ifdef FORCEFREE + DBPutQuadvar1(file, "uint_ff","mesh1", uintff, + dimensions, ndim, NULL, 0, + DB_DOUBLE, DB_ZONECENT, optList); + + DBPutQuadvar1(file, "u0_perp","mesh1", u0_perp, + dimensions, ndim, NULL, 0, + DB_DOUBLE, DB_ZONECENT, optList); + DBPutQuadvar1(file, "u1_perp","mesh1", u1_perp, + dimensions, ndim, NULL, 0, + DB_DOUBLE, DB_ZONECENT, optList); + DBPutQuadvar1(file, "u2_perp","mesh1", u2_perp, + dimensions, ndim, NULL, 0, + DB_DOUBLE, DB_ZONECENT, optList); + DBPutQuadvar1(file, "u3_perp","mesh1", u3_perp, + dimensions, ndim, NULL, 0, + DB_DOUBLE, DB_ZONECENT, optList); +#endif + DBPutQuadvar1(file, "lorentz_perp","mesh1", lorentz_perp, + dimensions, ndim, NULL, 0, + DB_DOUBLE, DB_ZONECENT, optList); + DBPutQuadvar1(file, "lorentz_par","mesh1", lorentz_par, + dimensions, ndim, NULL, 0, + DB_DOUBLE, DB_ZONECENT, optList); + + DBPutQuadvar1(file, "muBe","mesh1", muBe, dimensions, ndim, NULL, 0, DB_DOUBLE, DB_ZONECENT, optList); @@ -1911,6 +1967,7 @@ int fprint_silofile(ldouble time, int num, char* folder, char* prefix) free(rho); free(uint); + free(uintff); free(entr); free(temp); free(Omega); @@ -1924,6 +1981,12 @@ int fprint_silofile(ldouble time, int num, char* folder, char* prefix) free(u2); free(u3); free(lorentz); + free(u0_perp); + free(u1_perp); + free(u2_perp); + free(u3_perp); + free(lorentz_perp); + free(lorentz_par); free(vx); free(vy); free(vz); diff --git a/u2p.c b/u2p.c index 8547806..63a7cb6 100644 --- a/u2p.c +++ b/u2p.c @@ -14,9 +14,6 @@ static FTYPE compute_specificentropy_wmrho0_idealgas(FTYPE rho0, FTYPE wmrho0,FT static FTYPE compute_dspecificSdwmrho0_wmrho0_idealgas(FTYPE rho0, FTYPE wmrho0,FTYPE gamma); static FTYPE compute_dspecificSdrho_wmrho0_idealgas(FTYPE rho0, FTYPE wmrho0, FTYPE gamma); static int f_u2p_entropy(ldouble Wp, ldouble* cons, ldouble *f, ldouble *df, ldouble *err,ldouble pgamma); - -static ldouble u2p_solver_ff_parallel(ldouble Wguess,ldouble* cons, int verbose); -static int f_u2p_parallel(ldouble W, ldouble* cons,ldouble *f,ldouble *df,ldouble *err); //********************************************************************** //calculates primitives in given cell basing on global array u[] @@ -161,7 +158,7 @@ u2p(ldouble *uu0, ldouble *pp, void *ggg, int corrected[3], int fixups[2], int t int verbose=1; int hdcorr=0; int radcor=0; - corrected[0]=corrected[1]=0; + corrected[0]=corrected[1]=corrected[2]=0; fixups[0]=fixups[1]=0; int u2pret,u2pentrret,ret; @@ -179,23 +176,25 @@ u2p(ldouble *uu0, ldouble *pp, void *ggg, int corrected[3], int fixups[2], int t ret=0; u2pret=-1; - //negative uu[RHO] = rho u^t + //negative uu[RHO] = rho u^t -> set to floor if(uu[RHO] * gdetu_inv < 0.) { int gix,giy,giz; mpi_local2globalidx(geom->ix,geom->iy,geom->iz,&gix,&giy,&giz); - if(verbose) printf("%4d > %4d %4d %4d > NEGUU > neg rhout - requesting fixup\n",PROCID,gix,giy,giz); - pp[RHO]=RHOFLOOR; //used when not fixing up + if(verbose) + printf("%4d > %4d %4d %4d > neg rhout - requesting fixup\n",PROCID,gix,giy,giz); + + pp[RHO]=RHOFLOOR; uu[RHO]=RHOFLOOR*gdetu; - ret=-2; // to request fixup in all cases + ret=-2; // to request fixup in all cases -#ifndef SWAPPAPC + #ifndef SWAPPAPC global_int_slot[GLOBALINTSLOT_NTOTALMHDFIXUPS]++; // count as fixup -#endif + #endif } - -#ifdef ENFORCEENTROPY - u2pret = -1; //skip hot energy-conserving inversion and go to entropy inversion + +#if defined(ENFORCEENTROPY) || defined(FORCEFREE) + u2pret = -1; //skip hot energy-conserving inversion and go straight to entropy inversion #else u2pret = u2p_solver(uu,pp,ggg,U2P_HOT,0); // invert using the hot energy equation #endif //ENFORCEENTROPY @@ -208,7 +207,8 @@ u2p(ldouble *uu0, ldouble *pp, void *ggg, int corrected[3], int fixups[2], int t if(verbose>2 ) { - printf("u2p_entr >>> %d %d <<< %d >>> %e > %e\n",geom->ix + TOI, geom->iy + TOJ,u2pret,u0,pp[UU]); + printf("u2p_entr %d > %d %d %d \n", + u2pret, geom->ix + TOI, geom->iy + TOJ,geom->iz + TOK); } //************************************ @@ -221,7 +221,8 @@ u2p(ldouble *uu0, ldouble *pp, void *ggg, int corrected[3], int fixups[2], int t if(verbose>1) { - printf("u2p_entr err No. %4d > %e %e %e > %e %e > %4d %4d %4d\n",u2pret,uu[RHO],uu[UU],uu[ENTR],pp[RHO],pp[UU],geom->ix,geom->iy,geom->iz); + printf("u2p_entr err %4d > %4d %4d %4d\n", + u2pret,geom->ix+TOI,geom->iy+TOJ,geom->iz+TOK); } } // if(u2pret<0) // second time -- entropy eqn @@ -231,9 +232,10 @@ u2p(ldouble *uu0, ldouble *pp, void *ggg, int corrected[3], int fixups[2], int t if(u2pret<0) // entropy equation also failed { //leaving primitives unchanged - if(verbose > 1 || 1) + if(verbose > 1 || 1) // TODO ANDREW { - printf("%4d > %4d %4d %4d > MHDU2PFAIL > u2p prim. unchanged > %d \n",PROCID,geom->ix+TOI,geom->iy+TOJ,geom->iz+TOK,u2pret); + printf("%4d > MHDU2PFAIL %d > %4d %4d %4d > u2p prim. unchanged \n", + PROCID,u2pret,geom->ix+TOI,geom->iy+TOJ,geom->iz+TOK); } ret=-3; for(u2pret=0;u2pretix+TOI,geom->iy+TOJ,geom->iz+TOK,geom->ix,geom->iy,geom->iz,pp[RHO],TI,TJ,TK); - pp[RHO]=RHOFLOOR; + if(verbose) + { + printf("hd_floors CASE 1 at %d %d %d (%e) \n", + geom->ix+TOI,geom->iy+TOJ,geom->iz+TOK,pp[RHO]); + } + + pp[RHO]=RHOFLOOR; - ret=-1; + ret=-1; } //********************************************************************** - //rho too small, + // rho too small, // floor scaling as r^-3/2 #ifdef RHOFLOOR_BH ldouble xxBL[4]; @@ -370,16 +376,21 @@ check_floors_mhd(ldouble *pp, int whichvel,void *ggg) pporg[iv]=pp[iv]; } - if(verbose ) printf("hd_floors BH CASE 1 at %d %d (%e)\n",geom->ix+TOI,geom->iy+TOJ,pp[RHO]); - pp[RHO]=rhofloor; + if(verbose) + { + printf("hd_floors BH CASE 1 at %d %d %d (%e)\n", + geom->ix+TOI,geom->iy+TOJ,geom->iz+TOK,pp[RHO]); + } - ret=-1; + pp[RHO]=rhofloor; + + ret=-1; } #endif //*********************************************************************** - //rho too small, BH-disk like: - //Here we use the initial atmosphere as the floor on both density and pressure + // rho too small, BH-disk like: + // Here we use the initial atmosphere as the floor on both density and pressure #ifdef RHOFLOOR_INIT ldouble xxBL[4], rout = 2.; @@ -399,7 +410,12 @@ check_floors_mhd(ldouble *pp, int whichvel,void *ggg) pporg[iv] = pp[iv]; } - if(verbose ) printf("hd_floors BH CASE INIT at %d %d (%e)\n",geom->ix+TOI, geom->iy+TOJ, pp[RHO]); + if(verbose) + { + printf("hd_floors INIT CASE 1 at %d %d %d (%e)\n", + geom->ix+TOI, geom->iy+TOJ, geom->iz+TOK,pp[RHO]); + } + pp[RHO] = rhofloor; pp[UU] = uintfloor; @@ -416,10 +432,15 @@ check_floors_mhd(ldouble *pp, int whichvel,void *ggg) pporg[iv]=pp[iv]; } - if(verbose) {printf("hd_floors CASE 2 at (%d,%d,%d | %d,%d,%d): %e %e | tijk: %d %d %d\n",geom->ix+TOI,geom->iy+TOJ,geom->iz+TOK,geom->ix,geom->iy,geom->iz,pp[RHO],pp[UU],TI,TJ,TK);}//getchar();} - pp[UU]=UURHORATIOMIN*pp[RHO]; //increasing uint + if(verbose) + { + printf("hd_floors CASE 2 at %d,%d,%d (%e %e) \n", + geom->ix+TOI,geom->iy+TOJ,geom->iz+TOK,pp[RHO],pp[UU]); + //getchar(); + } - + pp[UU] = UURHORATIOMIN*pp[RHO]; //increasing uint + ret=-1; } @@ -432,10 +453,17 @@ check_floors_mhd(ldouble *pp, int whichvel,void *ggg) pporg[iv]=pp[iv]; } - pp[UU]=UURHORATIOMAX*pp[RHO]; //decreasing uint + if(verbose) + { + printf("hd_floors CASE 3 at %d,%d,%d (%e %e)\n", + geom->ix+TOI,geom->iy+TOJ,geom->iz+TOK,pp[RHO],pp[UU]); + //getchar(); + } + + pp[UU] = UURHORATIOMAX*pp[RHO]; //decreasing uint ret=-1; - if(verbose ) printf("hd_floors CASE 3 at (%d,%d,%d): %e %e\n",geom->ix+TOI,geom->iy+TOJ,geom->iz,pp[RHO],pp[UU]); + } //********************************************************************** @@ -451,6 +479,7 @@ check_floors_mhd(ldouble *pp, int whichvel,void *ggg) for(iv=1;iv<4;iv++) ucon[iv]=pp[1+iv]; + calc_ucon_ucov_from_prims(pp, geom, ucon, ucov); calc_bcon_bcov_bsq_from_4vel(pp, ucon, ucov, geom, bcon, bcov, &bsq); magpre = 0.5 * bsq; @@ -459,7 +488,11 @@ check_floors_mhd(ldouble *pp, int whichvel,void *ggg) // check density vs bsq if(bsq>B2RHORATIOMAX*pp[RHO]) { - if(verbose) printf("mag_floors CASE 1 at (%d,%d,%d): %e %e\n",geom->ix+TOI,geom->iy+TOJ,geom->iz,pp[RHO],bsq); + if(verbose) + { + printf("mag_floors CASE 1 at %d,%d,%d (%e %e)\n", + geom->ix+TOI,geom->iy+TOJ,geom->iz,pp[RHO],bsq); + } f=bsq/(B2RHORATIOMAX*pp[RHO]); //increase factor, f>1 rhofloored=1; } @@ -467,16 +500,19 @@ check_floors_mhd(ldouble *pp, int whichvel,void *ggg) // check ugas vs bsq if(bsq>B2UURATIOMAX*pp[UU]) { - if(verbose) printf("mag_floors CASE 2 at (%d,%d,%d): %e %e\n",geom->ix+TOI,geom->iy+TOJ,geom->iz,pp[UU],bsq); + if(verbose) + { + printf("mag_floors CASE 2 at (%d,%d,%d): %e %e\n", + geom->ix+TOI,geom->iy+TOJ,geom->iz,pp[UU],bsq); + } fuu=bsq/(B2UURATIOMAX*pp[UU]); //increase factor, fuu>1 uufloored=1; - //pp[UU]*=bsq/(B2UURATIOMAX*pp[UU]); } // floors were activated, apply them if(rhofloored==1 || uufloored==1) { - // backups + // save backups for(iv=0;iv 1e-2 ) - { - printf("(%d,%d,%d): %g %g QdotBold = %g, QdotBnew = %g, diff = %g\n", - geom->ix,geom->iy,geom->iz,f,fuu, - QdotB, QdotBnew, reldiff); - } -#elif(B2RHOFLOORFRAME==ZAMOFRAME) //new mass in ZAMO +#elif(B2RHOFLOORFRAME==ZAMOFRAME) //new mass in ZAMO frame ldouble dpp[NV],duu[NV]; ldouble drho=pp[RHO]*(f-1.); @@ -574,9 +595,9 @@ check_floors_mhd(ldouble *pp, int whichvel,void *ggg) conv_vels_ut(etacon,etarel,VEL4,VELPRIM,gg,GG); dpp[RHO]=drho; - //ANDREW this is a change, used to keep UU fixed here + //ANDREW this is a change,we used to keep UU fixed here dpp[UU]=dugas; - //dpp[UU]=0.; + //dpp[UU]=0.; dpp[VX] = etarel[1]; dpp[VY] = etarel[2]; dpp[VZ] = etarel[3]; @@ -593,145 +614,143 @@ check_floors_mhd(ldouble *pp, int whichvel,void *ggg) // find new prims after adding delta conserved in ZAMO int rettemp=0; - rettemp=u2p_solver(uu,pp,geom,U2P_HOT,0); + rettemp = u2p_solver(uu,pp,geom,U2P_HOT,0); if(rettemp<0) - rettemp=u2p_solver(uu,pp,geom,U2P_ENTROPY,0); + rettemp = u2p_solver(uu,pp,geom,U2P_ENTROPY,0); if(rettemp<0) // new mass in ZAMO frame failed { - printf("u2p failed after imposing bsq over rho floors at %d %d %d\n",geom->ix+TOI,geom->iy+TOJ,geom->iz+TOK); + printf("u2p failed after imposing bsq over rho floors at %d %d %d\n", + geom->ix+TOI,geom->iy+TOJ,geom->iz+TOK); -#ifdef B2RHOFLOOR_BACKUP_FFFRAME - // Backup bsq/rho floor -- if zamo frame fails, do fluid frame instead of crashing + #ifdef B2RHOFLOOR_BACKUP_FFFRAME + // if zamo frame fails, do fluid frame instead of crashing for(iv=0;ivix,geom->iy,geom->iz); + //keep energy density in ions and electrons fixed after applying the floors + ldouble Tg,Te,Ti; + Tg=calc_PEQ_Teifrompp(pporg,&Te,&Ti,geom->ix,geom->iy,geom->iz); - ldouble ne=calc_thermal_ne(pporg); //thermal only - ldouble ni=pporg[RHO]/MU_I/M_PROTON; + ldouble ne=calc_thermal_ne(pporg); //thermal only + ldouble ni=pporg[RHO]/MU_I/M_PROTON; - ldouble pe=K_BOLTZ*ne*Te; - ldouble pi=K_BOLTZ*ni*Ti; + ldouble pe=K_BOLTZ*ne*Te; + ldouble pi=K_BOLTZ*ni*Ti; - ldouble gammae=GAMMAE; - ldouble gammai=GAMMAI; - #ifdef CONSISTENTGAMMA - #ifndef FIXEDGAMMASPECIES - gammae=calc_gammaintfromtemp(Te,ELECTRONS); - gammai=calc_gammaintfromtemp(Ti,IONS); - #endif - #endif - - ldouble ue=pe/(gammae-1.); - ldouble ui=pi/(gammai-1.); - - //calculate new entropy using pp[] - ldouble mass, gamma, Tenew,Tinew, Senew, Sinew; - ldouble n, pnew; - - //electrons - mass = M_ELECTR; - gamma = GAMMAE; - n = calc_thermal_ne(pp); - #ifdef CONSISTENTGAMMA - Tenew=solve_Teifromnmu(n, mass, ue,ELECTRONS); //solves in parallel for gamma and temperature - theta=K_BOLTZ*Tenew/mass; - #ifndef FIXEDGAMMASPECIES - gamma=calc_gammaintfromtheta(theta); //the same gamma as just solved - #endif - #endif - pnew=(ue)*(gamma-1.); - Tenew=pnew/K_BOLTZ/n; - Senew=calc_SefromrhoT(n*MU_E*M_PROTON,Tenew,ELECTRONS); - pp[ENTRE]=Senew; - - //ions - mass = M_PROTON; - gamma = GAMMAI; - n = pp[RHO]/MU_I/M_PROTON;//calc_thermal_ne(pp); - #ifdef CONSISTENTGAMMA - Tinew=solve_Teifromnmu(n, mass, ui,IONS); //solves in parallel for gamma and temperature - theta=K_BOLTZ*Tinew/mass; - #ifndef FIXEDGAMMASPECIES - gamma=calc_gammaintfromtheta(theta); //the same gamma as just solved - #endif - #endif - pnew=(ui)*(gamma-1.); - Tinew=pnew/K_BOLTZ/n; - Sinew=calc_SefromrhoT(pp[RHO],Tinew,IONS); - pp[ENTRI]=Sinew; + ldouble gammae=GAMMAE; + ldouble gammai=GAMMAI; + #ifdef CONSISTENTGAMMA + #ifndef FIXEDGAMMASPECIES + gammae=calc_gammaintfromtemp(Te,ELECTRONS); + gammai=calc_gammaintfromtemp(Ti,IONS); + #endif + #endif + + ldouble ue=pe/(gammae-1.); + ldouble ui=pi/(gammai-1.); + + //calculate new entropy using updated pp[] + ldouble mass, gamma, theta; + ldouble Tenew,Tinew, Senew, Sinew; + ldouble n, pnew; + + //electrons + mass = M_ELECTR; + gamma = GAMMAE; + n = calc_thermal_ne(pp); + #ifdef CONSISTENTGAMMA + Tenew=solve_Teifromnmu(n, mass, ue,ELECTRONS); //solves in parallel for gamma and temperature + theta=K_BOLTZ*Tenew/mass; + #ifndef FIXEDGAMMASPECIES + gamma=calc_gammaintfromtheta(theta); //the same gamma as just solved + #endif + #endif + pnew=(ue)*(gamma-1.); + Tenew=pnew/K_BOLTZ/n; + Senew=calc_SefromrhoT(n*MU_E*M_PROTON,Tenew,ELECTRONS); + pp[ENTRE]=Senew; + + //ions + mass = M_PROTON; + gamma = GAMMAI; + n = pp[RHO]/MU_I/M_PROTON;//calc_thermal_ne(pp); + #ifdef CONSISTENTGAMMA + Tinew=solve_Teifromnmu(n, mass, ui,IONS); //solves in parallel for gamma and temperature + theta=K_BOLTZ*Tinew/mass; + #ifndef FIXEDGAMMASPECIES + gamma=calc_gammaintfromtheta(theta); //the same gamma as just solved + #endif + #endif + pnew=(ui)*(gamma-1.); + Tinew=pnew/K_BOLTZ/n; + Sinew=calc_SefromrhoT(pp[RHO],Tinew,IONS); + pp[ENTRI]=Sinew; #endif //EVOLVEELECTRONS - ret=-1; - } //if(bsq>B2RHORATIOMAX*pp[RHO]) or (bsq>B2UURATIOMAX*pp[UU]) + ret=-1; + + } // if(rhofloored==1 || uufloored==1) #endif //MAGNFIELD //********************************************************************** //too fast - //TODO: implement checks for other VELPRIM + //TODO: implement GAMMAMAXHD for other VELPRIM if(VELPRIM==VELR) { - ldouble qsq=0.; - int i,j; - for(i=1;i<4;i++) - for(j=1;j<4;j++) - qsq+=pp[UU+i]*pp[UU+j]*gg[i][j]; - ldouble gamma2=1.+qsq; - if(gamma2>GAMMAMAXHD*GAMMAMAXHD) + ldouble qsq=0.; + int i,j; + for(i=1;i<4;i++) + for(j=1;j<4;j++) + qsq+=pp[UU+i]*pp[UU+j]*gg[i][j]; + ldouble gamma2=1.+qsq; + + if(gamma2>GAMMAMAXHD*GAMMAMAXHD) + { + ldouble qsqmax=GAMMAMAXHD*GAMMAMAXHD-1.; + ldouble A=sqrt(qsqmax/qsq); + for(j=1;j<4;j++) + pp[UU+j]*=A; + + if(verbose) { - ldouble qsqmax=GAMMAMAXHD*GAMMAMAXHD-1.; - ldouble A=sqrt(qsqmax/qsq); - for(j=1;j<4;j++) - pp[UU+j]*=A; - ret=-1; - if(verbose ) - { - printf("hd_floors CASE 4 at (%d,%d,%d): %e",geom->ix+TOI,geom->iy+TOJ,geom->iz,sqrt(gamma2)); - qsq=0.; - for(i=1;i<4;i++) - for(j=1;j<4;j++) - qsq+=pp[UU+i]*pp[UU+j]*gg[i][j]; - gamma2=1.+qsq; - printf(" -> %e\n",sqrt(gamma2)); - } + printf("hd_floors CASE 4 at %d,%d,%d (%e)",geom->ix+TOI,geom->iy+TOJ,geom->iz,sqrt(gamma2)); } + } + + ret = -1; } - - //********************************************************************** //Species temperature floors/ceilings + //TODO ANDREW will not keep them consistent ue + ui = ugas! #ifdef EVOLVEELECTRONS ldouble mue,mui; mui=MU_I; @@ -780,7 +799,7 @@ check_floors_mhd(ldouble *pp, int whichvel,void *ggg) ldouble Timinimal=TEMPIMINIMALFRACTION*Tgas; if(Tilocix,geom->iy,geom->iz); - + #ifdef FORCEFREE + //ldouble uutmp[NV]; // ANDREW TODO why initial crash??? + pp[UUFF] = pp[UU]; + #endif + } +#endif //SKIPALLFLOORS + return ret; } //********************************************************************** +// solver wrapper //********************************************************************** -//********************************************************************** -//routines for various types of inversions -//********************************************************************** -//********************************************************************** -//********************************************************************** - -//******************************************** -//Harm u2p_hot -//******************************************** - -static FTYPE -dpdWp_calc_vsq(FTYPE Wp, FTYPE D, FTYPE vsq, FTYPE gamma) -{ - FTYPE W=Wp+D; - return( (gamma - 1.) * (1. - vsq) / gamma ) ; -} - -// 1 / (d(u+p)/dp) -static FTYPE -compute_idwmrho0dp(FTYPE wmrho0, FTYPE gamma) -{ - return((gamma-1.)/gamma); -} - - -// 1 / (drho0/dp) holding wmrho0 fixed -static FTYPE -compute_idrho0dp(FTYPE wmrho0) -{ - return(0.0); -} +//ANDREW TODO do we need to use function pointers? -static int -f_u2p_hot(ldouble Wp, ldouble* cons,ldouble *f,ldouble *df,ldouble *err,ldouble pgamma) +int +u2p_solver(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose) { + struct geometry *geom + = (struct geometry *) ggg; - ldouble Qn=cons[0]; - ldouble Qt2=cons[1]; - ldouble D=cons[2]; - ldouble QdotBsq=cons[3]; - ldouble Bsq=cons[4]; - ldouble Qdotnp=cons[6]; - - ldouble W=Wp+D; - - FTYPE W3,X3,Ssq,Wsq,X,X2,Xsq; - FTYPE Qtsq = Qt2; - X = Bsq + W; - Wsq = W*W; - W3 = Wsq*W ; - X2 = X*X; - Xsq = X2; - X3 = X2*X; - - ldouble v2=( Wsq * Qtsq + QdotBsq * (Bsq + 2.*W)) / (Wsq*Xsq); - ldouble gamma2 = 1./(1.-v2); - ldouble gamma = sqrt(gamma2); - ldouble rho0 = D/gamma; - ldouble wmrho0 = Wp/gamma2 - D*v2/(1.+gamma); - ldouble u = wmrho0 / pgamma; - ldouble p = (pgamma-1)*u; - - //original: -#if (U2P_EQS==U2P_EQS_NOBLE) - *f = Qn + W - p + 0.5*Bsq*(1.+v2) - QdotBsq/2./Wsq; - *err = fabs(*f) / (fabs(Qn) + fabs(W) + fabs(p) + fabs(0.5*Bsq*(1.+v2)) + fabs(QdotBsq/2./Wsq)); -#endif - - //JONS: -#if (U2P_EQS==U2P_EQS_JON) - *f = Qdotnp + Wp - p + 0.5*Bsq + (Bsq*Qtsq - QdotBsq)/X2; - *err = fabs(*f) / (fabs(Qdotnp) + fabs(Wp) + fabs(p) + fabs(0.5*Bsq) + fabs((Bsq*Qtsq - QdotBsq)/X2)); +#ifdef NONRELMHD + return u2p_solver_nonrel(uu,pp,ggg,Etype,verbose); #endif - // dp/dWp = dp/dW + dP/dv^2 dv^2/dW - ldouble dvsq=(-2.0/X3 * ( Qtsq + QdotBsq * (3.0*W*X + Bsq*Bsq)/W3)); - ldouble dp1 = dpdWp_calc_vsq(Wp, D, v2 ,pgamma); // vsq can be unphysical - - ldouble idwmrho0dp=compute_idwmrho0dp(wmrho0,pgamma); - ldouble dwmrho0dvsq = (D*(gamma*0.5-1.0) - Wp); - - ldouble drho0dvsq = -D*gamma*0.5; // because \rho = D/\gamma - ldouble idrho0dp = compute_idrho0dp(wmrho0); - - ldouble dp2 = drho0dvsq *idrho0dp + dwmrho0dvsq *idwmrho0dp; - - ldouble dpdW = dp1 + dp2*dvsq; // dp/dW = dp/dWp - - //original: - #if (U2P_EQS==U2P_EQS_NOBLE) - *df=1.-dpdW + QdotBsq/(Wsq*W) + 0.5*Bsq*dvsq; - #endif - - //JONS: - #if (U2P_EQS==U2P_EQS_JON) - *df=1. -dpdW + (Bsq*Qtsq - QdotBsq)/X3*(-2.0); + int do_u2p_ff=0; +#ifdef FORCEFREE + #ifndef HYBRID_FORCEFREE + do_u2p_ff=1; + #else + + //ANDREW TODO what should be the criterion here? local sigma from guess pp? + ldouble ucon0[4],ucov0[4],bcon0[4],bcov0[4]; + ldouble bsq0,sigma0; + calc_ucon_ucov_from_prims(pp, geom, ucon0, ucov0); + calc_bcon_bcov_bsq_from_4vel(pp, ucon0, ucov0, geom, bcon0, bcov0, &bsq0); + sigma0 = bsq0/pp[RHO]; + if(sigma0 > HYBRID_FORCEFREE_SIGMACUT) + do_u2p_ff = 1; #endif +#endif //FORCEFREE - return 0; -} - -//******************************************** -//Harm u2p_entropy -//******************************************** - -// p(rho0, w-rho0 = u+p) -static FTYPE -pressure_wmrho0_idealgas(FTYPE rho0, FTYPE wmrho0,FTYPE gamma) -{ - ldouble igammar = (gamma-1.)/gamma; - return(igammar*wmrho0) ; -} - -// local aux function -static FTYPE -compute_inside_entropy_wmrho0_idealgas(FTYPE rho0, FTYPE wmrho0,FTYPE gamma) -{ - FTYPE pressure,indexn,insideentropy; - - pressure=pressure_wmrho0_idealgas(rho0,wmrho0,gamma); - indexn=1.0/(gamma-1.); - insideentropy=pow(pressure,indexn)/pow(rho0,indexn+1.0); - - return(insideentropy); -} - - -// specific entropy as function of rho0 and internal energy (u) -// Ss(rho0,\chi=u+p) -// specific entropy = \ln( p^n/\rho^{n+1} ) -static FTYPE -compute_specificentropy_wmrho0_idealgas(FTYPE rho0, FTYPE wmrho0,FTYPE gamma) -{ - FTYPE insideentropy,specificentropy; - - insideentropy=compute_inside_entropy_wmrho0_idealgas(rho0, wmrho0,gamma); - specificentropy=log(insideentropy); - - return(specificentropy); - -} - -// used for utoprim_jon when doing entropy evolution -// dSspecific/d\chi -static FTYPE -compute_dspecificSdwmrho0_wmrho0_idealgas(FTYPE rho0, FTYPE wmrho0,FTYPE gamma) -{ - FTYPE dSdchi; + // use ff solver if not hybrid or if high sigma + if(do_u2p_ff) + { + return u2p_solver_ff(uu,pp,ggg,U2P_ENTROPY,verbose); + } + else + { + + //int (*solver)(ldouble*,ldouble*,void*,int,int); + + #if (U2P_SOLVER==U2P_SOLVER_WP) + //solver = & u2p_solver_Wp; + return u2p_solver_Wp(uu,pp,ggg,Etype,verbose); + #endif + + #if (U2P_SOLVER==U2P_SOLVER_W) // this is the default + //solver = & u2p_solver_W; + return u2p_solver_W(uu,pp,ggg,Etype,verbose); + #endif + } - dSdchi = 1.0/((gamma-1.)*wmrho0); - // Again, GAMMA->1 means dSdchi->\infty unless \chi->0 or rho0->0 + return -200; + //return (*solver)(uu,pp,ggg,Etype,verbose); +} - return(dSdchi); -} +//********************************************************************** +//Newton-Raphson solver +//iterates W, not Wp +//Etype == 0 -> hot inversion (uses D,Ttt,Tti) +//Etype == 1 -> entropy inversion (uses D,S,Tti) +//********************************************************************** -// dSspecific/drho0 -static FTYPE -compute_dspecificSdrho_wmrho0_idealgas(FTYPE rho0, FTYPE wmrho0, FTYPE gamma) +int +u2p_solver_W(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose) { - FTYPE dSdrho; + int i,j,k; + ldouble rho,uint,w,W,alpha,D,Sc; + ldouble ucon[4],ucov[4],utcon[4],utcov[4],ncov[4],ncon[4]; + ldouble Qcon[4],Qcov[4],Qconp[4],Qcovp[4],jmunu[4][4],Qtcon[4],Qtcov[4],Qt2,Qn; + ldouble QdotB,QdotBsq,Bcon[4],Bcov[4],Bsq; - dSdrho=gamma/((1.0-gamma)*rho0); - - return(dSdrho); -} - -static int -f_u2p_entropy(ldouble Wp, ldouble* cons, ldouble *f, ldouble *df, ldouble *err,ldouble pgamma) -{ - ldouble Qn=cons[0]; - ldouble Qt2=cons[1]; - ldouble D=cons[2]; - ldouble QdotBsq=cons[3]; - ldouble Bsq=cons[4]; - ldouble Sc=cons[5]; - - ldouble W=Wp+D; - - FTYPE W3,X3,Ssq,Wsq,X,X2,Xsq; - FTYPE Qtsq = Qt2; - X = Bsq + W; - Wsq = W*W; - W3 = Wsq*W ; - X2 = X*X; - Xsq = X2; - X3 = X2*X; - - ldouble v2=( Wsq * Qtsq + QdotBsq * (Bsq + 2.*W)) / (Wsq*Xsq); - ldouble gamma2 = 1./(1.-v2); - ldouble gamma = sqrt(gamma2); - ldouble rho0 = D/gamma; - ldouble wmrho0 = Wp/gamma2 - D*v2/(1.+gamma); - ldouble u = wmrho0 / pgamma; - ldouble p = (pgamma-1)*u; - - ldouble Ssofchi=compute_specificentropy_wmrho0_idealgas(rho0,wmrho0,pgamma); - - *f= -Sc + D*Ssofchi; - - *err = fabs(*f) / (fabs(Sc) + fabs(D*Ssofchi)); - - FTYPE dSsdW,dSsdvsq,dSsdWp,dScprimedWp,dSsdrho,dSsdchi; - FTYPE dvsq,dwmrho0dW,drho0dW; - FTYPE dwmrho0dvsq,drho0dvsq; - - dSsdrho=compute_dspecificSdrho_wmrho0_idealgas(rho0,wmrho0,pgamma); - dSsdchi=compute_dspecificSdwmrho0_wmrho0_idealgas(rho0,wmrho0,pgamma); - - dwmrho0dW = 1.0/gamma; // holding utsq fixed - drho0dW = 0.0; // because \rho=D/\gamma and holding utsq fixed - dwmrho0dvsq = (D*(gamma*0.5-1.0) - Wp); // holding Wp fixed - drho0dvsq = -D*gamma*0.5; // because \rho=D/\gamma and holding Wp fixed - - dvsq=(-2.0/X3 * ( Qtsq + QdotBsq * (3.0*W*X + Bsq*Bsq)/W3)); - - dSsdW = drho0dW *dSsdrho + dwmrho0dW *dSsdchi; // dSs/dW' holding utsq fixed - dSsdvsq = drho0dvsq *dSsdrho + dwmrho0dvsq *dSsdchi; - dSsdWp = dSsdW + dSsdvsq*dvsq; // dSs/dW = dSs/dWp [total derivative] - - dScprimedWp = D*dSsdWp; - - *df = dScprimedWp; - - return 0; - -} - -//********************************************************************** -// solver wrapper -//********************************************************************** -//ANDREW TODO do we need function pointers? - -int -u2p_solver(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose) -{ - -#ifdef NONRELMHD - return u2p_solver_nonrel(uu,pp,ggg,Etype,verbose); -#endif - -#ifdef FORCEFREE - return u2p_solver_ff(uu,pp,ggg,1); -#endif - - int (*solver)(ldouble*,ldouble*,void*,int,int); - struct geometry *geom - = (struct geometry *) ggg; - -#if (U2P_SOLVER==U2P_SOLVER_WP) - solver = & u2p_solver_Wp; -#endif - -#if (U2P_SOLVER==U2P_SOLVER_W) // this is the default - solver = & u2p_solver_W; -#endif - - return (*solver)(uu,pp,ggg,Etype,verbose); -} - -//********************************************************************** -//non-relativistic, analytical, solver -//********************************************************************** - -int -u2p_solver_nonrel(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose) -{ + /*******************************************************/ //prepare geometry struct geometry *geom - = (struct geometry *) ggg; - - ldouble (*gg)[5], (*GG)[5], gdet, gdetu, gdetu_inv; - gg=geom->gg; GG=geom->GG; - gdet=geom->gdet;gdetu=gdet; -#if (GDETIN==0) //gdet out of derivatives - gdetu=1.; -#endif - gdetu_inv = 1. / gdetu; - - //density - ldouble rho=uu[RHO] * gdetu_inv; - pp[RHO]=rho; - - //velocities u_i - ldouble ucov[4],ucon[4],vcov[4]; - ucov[0]=-1.; - ucov[1]=uu[VX] * gdetu_inv / rho; - ucov[2]=uu[VY] * gdetu_inv / rho; - ucov[3]=uu[VZ] * gdetu_inv / rho; - - indices_12(ucov,ucon,GG); - - ucon[0]=1.; - - ldouble v2=dot3nr(ucon,ucov); - - pp[VX]=ucon[1]; - pp[VY]=ucon[2]; - pp[VZ]=ucon[3]; - - ldouble bsq=0.; - -#ifdef MAGNFIELD - - ldouble bcon[4],bcov[4],Bcon[4]; - - Bcon[0]=0.; - Bcon[1]=uu[B1] * gdetu_inv; - Bcon[2]=uu[B2] * gdetu_inv ; - Bcon[3]=uu[B3] * gdetu_inv ; - - pp[B1]=Bcon[1]; - pp[B2]=Bcon[2]; - pp[B3]=Bcon[3]; - - calc_bcon_bcov_bsq_from_4vel(pp, ucon, ucov, geom, bcon, bcov, &bsq); - -#endif - + = (struct geometry *) ggg; - ldouble uint; - if(Etype==U2P_HOT) - { - uint = -uu[UU] * gdetu_inv - bsq/2. - rho*v2/2.; - - if(uint %e %e %e %e %e\n",geom->ix+TOI,geom->iy+TOI,uint,uu[UU]/gdetu,bsq/2.,rho*v2/2.,rho); - //printf("%e %e\n",uint,rho); - //if(geom->ix>50) getch(); - return -1; - } - - pp[UU]=uint; - } - else if(Etype==U2P_ENTROPY) - { - ldouble S=uu[ENTR] * gdetu_inv ; - uint= calc_ufromS(S,rho,geom->ix,geom->iy,geom->iz); - pp[UU]=uint; - } - - //pure entropy evolution - updated only in the end of RK2 - pp[ENTR]= uu[ENTR] * gdetu_inv; - - #ifdef EVOLVEELECTRONS - ldouble Se=uu[ENTRE] * gdetu_inv ; - pp[ENTRE]=Se; - - ldouble Si=uu[ENTRI] * gdetu_inv ; - pp[ENTRI]=Si; - #endif - - #ifdef RELELECTRONS - int ib; - for(ib=0;ib hot inversion (uses D,Ttt,Tti) -//Etype == 1 -> entropy inversion (uses D,S,Tti) -//********************************************************************** - -int -u2p_solver_Wp(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose) -{ - int i,j,k; - ldouble rho,uint,w,W,Wp,Wpprev,alpha,D,Sc,alphasq,betasqoalphasq; - ldouble ucon[4],ucov[4],utcon[4],utcov[4],ncov[4],ncon[4]; - ldouble Qcon[4],Qcov[4],Qconp[4],Qcovp[4],jmunu[4][4],Qtcon[4],Qtcov[4],Qt2,Qn,Qdotnp; - ldouble QdotB,QdotBsq,Bcon[4],Bcov[4],Bsq; - - /****************************/ - //prepare geometry - struct geometry *geom - = (struct geometry *) ggg; - ldouble pgamma=GAMMA; - #ifdef CONSISTENTGAMMA +#ifdef CONSISTENTGAMMA pgamma=pick_gammagas(geom->ix,geom->iy,geom->iz); - #endif +#endif ldouble pgammam1=pgamma-1.; - + ldouble (*gg)[5], (*GG)[5], gdet, gdetu, gdetu_inv; - gg=geom->gg; GG=geom->GG; - gdet=geom->gdet;gdetu=gdet; + gg=geom->gg; GG=geom->GG; + gdet=geom->gdet; gdetu=gdet; #if (GDETIN==0) //gdet out of derivatives gdetu=1.; #endif gdetu_inv = 1. / gdetu; - + /****************************/ //equations choice int (*f_u2p)(ldouble,ldouble*,ldouble*,ldouble*,ldouble*,ldouble); - if(Etype==U2P_HOT) - f_u2p=&f_u2p_hot; - if(Etype==U2P_ENTROPY) - f_u2p=&f_u2p_entropy; + if(Etype==U2P_HOT) + f_u2p=&f_u2p_hot; + if(Etype==U2P_ENTROPY) + f_u2p=&f_u2p_entropy; /****************************/ - + if(verbose>1) { printf("********************\n"); print_conserved(uu); print_primitives(pp); } - + //conserved quantities etc - //alpha - alpha=geom->alpha; - alphasq=alpha*alpha; - + //alpha=sqrt(-1./GG[0][0]); + alpha = geom->alpha; + //D D = uu[RHO] * gdetu_inv * alpha; // gdetu rho ut - - //conserved entropy "gdet S u^t" + + //conserved entropy "S u^t" Sc = uu[ENTR] * gdetu_inv * alpha; - - //Qp_mu=alpha T^t_mu - Qcovp[0]=uu[UU] * gdetu_inv * alpha; - Qcovp[1]=uu[VX] * gdetu_inv * alpha; - Qcovp[2]=uu[VY] * gdetu_inv * alpha; - Qcovp[3]=uu[VZ] * gdetu_inv * alpha; - + + //Q_mu=alpha T^t_mu + Qcov[0] = (uu[UU] * gdetu_inv - uu[RHO] * gdetu_inv) * alpha; + Qcov[1] = uu[VX] * gdetu_inv * alpha; + Qcov[2] = uu[VY] * gdetu_inv * alpha; + Qcov[3] = uu[VZ] * gdetu_inv * alpha; + + //Qp_mu=alpha T^t_mu + Qcovp[0] = uu[UU] * gdetu_inv *alpha; + Qcovp[1] = uu[VX] * gdetu_inv *alpha; + Qcovp[2] = uu[VY] * gdetu_inv *alpha; + Qcovp[3] = uu[VZ] * gdetu_inv *alpha; + //Qp^mu indices_12(Qcovp,Qconp,GG); - - //Q_mu=alpha (T^t_mu - rho u^t delta(t,mu)) - avoid this one - Qcov[0]=(uu[UU] * gdetu_inv - uu[RHO] * gdetu_inv ) * alpha; - Qcov[1]=uu[VX] * gdetu_inv * alpha; - Qcov[2]=uu[VY] * gdetu_inv * alpha; - Qcov[3]=uu[VZ] * gdetu_inv * alpha; - + //Q^mu indices_12(Qcov,Qcon,GG); - + #ifdef MAGNFIELD //curly B^mu Bcon[0]=0.; Bcon[1]=uu[B1] * gdetu_inv * alpha; Bcon[2]=uu[B2] * gdetu_inv * alpha; Bcon[3]=uu[B3] * gdetu_inv * alpha; - + //B_mu indices_21(Bcon,Bcov,gg); - + Bsq = dot(Bcon,Bcov); - QdotB = dot(Qcov,Bcon); - QdotBsq = QdotB*QdotB; - + #else Bsq=QdotB=QdotBsq=0.; Bcon[0]=Bcon[1]=Bcon[2]=Bcon[3]=0.; -#endif +#endif // MAGNFIELD + //normal observer velocity //n_mu = (-alpha, 0, 0, 0) ncov[0]=-alpha; ncov[1]=ncov[2]=ncov[3]=0.; //n^mu indices_12(ncov,ncon,GG); - + //Q_mu n^mu = Q^mu n_mu = -alpha*Q^t - Qn=Qcon[0] * ncov[0]; - - //\beta^i \beta_i / \alpha^2 = g^{ti} g_{ti} - betasqoalphasq=gg[0][1]*GG[0][1] + gg[0][2]*GG[0][2] + gg[0][3]*GG[0][3]; + Qn = Qcon[0] * ncov[0]; - //Qdotnp=-E'=-E+D - ldouble Dfactor = (-geom->gttpert + alphasq*betasqoalphasq)/(alphasq+alpha); - Qdotnp = Qconp[0]*ncov[0] + D*(Dfactor) ; // -Qdotn-W = -Qdotnp-Wp - - //j^mu_nu=delta^mu_nu +n^mu n_nu + //j^mu_nu = delta^mu_nu +n^mu n_nu for(i=0;i<4;i++) for(j=0;j<4;j++) jmunu[i][j] = delta(i,j) + ncon[i]*ncov[j]; - + //Qtilda^nu = j^nu_mu Q^mu for(i=0;i<4;i++) { @@ -1350,15 +1057,22 @@ u2p_solver_Wp(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose) for(j=0;j<4;j++) Qtcon[i]+=jmunu[i][j]*Qcon[j]; } - + //Qtilda_nu indices_21(Qtcon,Qtcov,gg); - + //Qt2=Qtilda^mu Qtilda_mu Qt2=dot(Qtcon,Qtcov); FTYPE Qtsq = Qt2; - - //initial guess for Wp = w gamma**2 based on primitives + + //\beta^i \beta_i / \alpha^2 = g^{ti} g_{ti} + ldouble betasqoalphasq=gg[0][1]*GG[0][1] + gg[0][2]*GG[0][2] + gg[0][3]*GG[0][3]; + ldouble alphasq=alpha*alpha; + //Qdotnp=-E'=-E+D + ldouble Dfactor = (-geom->gttpert + alphasq*betasqoalphasq)/(alphasq+alpha); + ldouble Qdotnp = Qconp[0]*ncov[0] + D*(Dfactor) ; // -Qdotn - W = -Qdotnp-Wp + + //initial guess for W = w gamma**2 based on current primitives rho=pp[RHO]; uint=pp[UU]; utcon[0]=0.; @@ -1370,203 +1084,175 @@ u2p_solver_Wp(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose) { conv_vels(utcon,utcon,VELPRIM,VELR,gg,GG); } - + ldouble qsq=0.; for(i=1;i<4;i++) for(j=1;j<4;j++) qsq+=utcon[i]*utcon[j]*gg[i][j]; ldouble gamma2=1.+qsq; ldouble gamma=sqrt(gamma2); - - //Wp - Wp=(pgamma*uint)*gamma2; - - ldouble Wpinit, Winit; - if(verbose>1) printf("initial Wp:%e\n",Wp); - Wpinit=Wp; - Winit=Wpinit+D; - - //test if does not provide reasonable gamma2 - // Make sure that W is large enough so that v^2 < 1 and w-rho > 0 : + + //W + W=(rho+pgamma*uint)*gamma2; + + if(verbose>1) printf("initial W:%e\n",W); + + // test if does not provide reasonable gamma2 + // Make sure that W is large enough so that v^2 < 1 : int i_increase = 0; ldouble f0,f1,dfdW,err; - ldouble CONV=U2PCONV; + ldouble CONV=U2PCONV; + ldouble EPS=1.e-4; + ldouble Wprev=W; ldouble cons[7]={Qn,Qt2,D,QdotBsq,Bsq,Sc,Qdotnp}; - int iter=0,fu2pret; do + { + f0=dfdW=0.; + + //if(Etype!=U2P_HOT) //entropy-like solvers require this additional check + //now invoked for all solvers: + (*f_u2p)(W-D,cons,&f0,&dfdW,&err,pgamma); + + if( ((( W*W*W * ( W + 2.*Bsq ) + - QdotBsq*(2.*W + Bsq) ) <= W*W*(Qtsq-Bsq*Bsq)) + || !isfinite(f0) || !isfinite(f0) + || !isfinite(dfdW) || !isfinite(dfdW)) + && (i_increase < 50)) { - W=Wp+D; - f0=dfdW=0.; - - FTYPE Wsq,Xsq,X; - X = Bsq + W; - Xsq = X*X; - Wsq = W*W; - - ldouble v2=( Wsq * Qtsq + QdotBsq * (Bsq + 2.*W)) / (Wsq*Xsq); - ldouble gamma2 = 1./(1.-v2); - ldouble gamma = sqrt(gamma2); - ldouble rho0 = D/gamma; - ldouble wmrho0 = Wp/gamma2 - D*v2/(1.+gamma); - - (*f_u2p)(Wp,cons,&f0,&dfdW,&err,pgamma); - - if((gamma2<0. || Wp<0. || wmrho0<0.|| !isfinite(f0) || !isfinite(dfdW)) && (i_increase < 50)) - { - if(verbose>0) printf("init Wp : %e - %e %e %e %e\n",Wp,v2,wmrho0,f0,dfdW); - Wp *= 2.; - i_increase++; - continue; - } - else - break; - } - while(1); - - if(i_increase>=50) - { - if(verbose>0) - {printf("failed to find initial W for Etype: %d\n",Etype); - printf("at %d %d\n",geom->ix+TOI,geom->iy+TOJ);} - return -150; - - print_NVvector(uu); - print_NVvector(pp); - //getchar(); - } - - //1d Newton solver - do - { - Wpprev=Wp; - iter++; - - fu2pret=(*f_u2p)(Wp,cons,&f0,&dfdW,&err,pgamma); - - if(verbose>1) printf("%d %e %e %e %e\n",iter,Wp,f0,dfdW,err); - - //convergence test - //if(err1) printf("sub (%d) :%d %e %e %e %e %e %e\n",idump,iter,Wpnew,f0tmp,dfdWtmp,v2,gamma2,wmrho0); - - //if((gamma2<0. || Wpnew<0. || wmrho0<0. || !isfinite(f0tmp) || !isfinite(dfdWtmp)) && (idump=itmaxdamp) - { - if(verbose>0) printf("damped unsuccessfuly\n"); - return -101; - } - - if(fabs(W)>BIG) - { - if(verbose>1) printf("W has gone out of bounds at %d,%d,%d\n",geom->ix+TOI,geom->iy+TOJ,geom->iz); - return -103; - } - - Wp=Wpnew; + if(verbose>0) printf("init W : %e -> %e (%e %e)\n",W,100.*W,f0,dfdW); + W *= 10.; + i_increase++; + continue; + } + else + break; + } + while(1); + + if(i_increase>=50) + { + return -150; + printf("failed to find initial W for Etype: %d\n",Etype); + printf("at %d %d\n",geom->ix+TOI,geom->iy+TOJ); + print_NVvector(uu); + print_NVvector(pp); + //getchar(); + } + + //1d Newton solver + int iter=0, fu2pret; + do + { + Wprev=W; + iter++; - //convergence test: - if(err1) printf("%d %e %e %e %e\n",iter,W,f0,dfdW,err); + + //convergence test + if(err1) printf("sub (%d) :%d %e %e %e %e\n",idump,iter,Wnew,f0tmp,dfdWtmp,errtmp); + if( ((( Wnew*Wnew*Wnew * ( Wnew + 2.*Bsq ) + - QdotBsq*(2.*Wnew + Bsq) ) <= Wnew*Wnew*(Qtsq-Bsq*Bsq)) + || !isfinite(f0tmp) || !isfinite(f0tmp) + || !isfinite(dfdWtmp) || !isfinite(dfdWtmp)) + && (idump<100)) + { + idump++; + dumpfac/=2.; + Wnew=W-dumpfac*f0/dfdW; + continue; + } + else + break; } - while(iter<50); - - if(iter>=50) + while(1); + + if(idump>=100) { - if(verbose>0) printf("iter exceeded in u2p_solver with Etype: %d\n",Etype); //getchar(); - return -102; + if(verbose>0) printf("damped unsuccessfuly\n"); + return -101; } - - if(!isfinite(Wp) || !isfinite(Wp)) {if(verbose) printf("nan/inf W in u2p_solver with Etype: %d\n",Etype); return -103;} - - if(verbose>1) + + W=Wnew; + + if(fabs(W)>BIG) { - fu2pret=(*f_u2p)(Wp,cons,&f0,&dfdW,&err,pgamma); - printf("end: %d %e %e %e %e\n",iter,Wp,f0,dfdW,err); + if(verbose>1) printf("W has gone out of bounds at %d,%d,%d\n",geom->ix+TOI,geom->iy+TOJ,geom->iz); + return -103; } - + + if(fabs((W-Wprev)/Wprev)=50) + { + if(verbose>0) printf("iter exceeded in u2p_solver with Etype: %d\n",Etype); //getchar(); + return -102; + } + + if(!isfinite(W)) + { + if(verbose) printf("nan/inf W in u2p_solver with Etype: %d\n",Etype); + return -103; + } + + if(verbose>1) + { + fu2pret=(*f_u2p)(W-D,cons,&f0,&dfdW,&err,pgamma); + printf("end: %d %e %e %e %e\n",iter,W,f0,dfdW,err); + } + //W found, let's calculate v2 and the rest - W=Wp+D; - ldouble Wsq,Xsq,v2,wmrho0,entr; - + ldouble Wsq,Xsq,v2,entr; + Wsq = W*W ; - Xsq = (Bsq + W) * (Bsq + W); + Xsq = (Bsq + W) * (Bsq + W); v2 = ( Wsq * Qtsq + QdotBsq * (Bsq + 2.*W)) / (Wsq*Xsq); - + gamma2=1./(1.-v2); gamma=sqrt(gamma2); rho=D/gamma; entr=Sc/gamma; - // w-\rho_0 = (u+p) = W'/\gamma^2 - D v^2/(1+\gamma) - wmrho0 = Wp/gamma2 - D*v2/(1.+gamma); - uint=1./pgamma*wmrho0; + uint=1./pgamma*(W/gamma2-rho); utcon[0]=0.; utcon[1]=gamma/(W+Bsq)*(Qtcon[1]+QdotB*Bcon[1]/W); utcon[2]=gamma/(W+Bsq)*(Qtcon[2]+QdotB*Bcon[2]/W); utcon[3]=gamma/(W+Bsq)*(Qtcon[3]+QdotB*Bcon[3]/W); + if(verbose>1) + printf("end2: %e %e %e %e %e %e\n",W,D,pgamma,gamma2,rho,uint); if(!isfinite(utcon[1])) - { - //print_4vector(utcon); - return -120; - } - - - if(uint<0. || gamma2<0. ||isnan(Wp) || !isfinite(Wp)) - { - if(verbose>0) printf("neg u in u2p_solver %e %e %e %e\n",rho,uint,gamma2,W);//getchar(); - return -104; - } - + { + return -120; + } + + if(uint<0. || gamma2<0. || isnan(W) || !isfinite(W)) + { + if(verbose>0) printf("neg u in u2p_solver %e %e %e %e\n",rho,uint,gamma2,W);//getchar(); + return -104; + } + //converting to VELPRIM if (VELR != VELPRIM) { @@ -1579,30 +1265,32 @@ u2p_solver_Wp(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose) pp[VX]=utcon[1]; pp[VY]=utcon[2]; pp[VZ]=utcon[3]; - - if(rho<0.) - { - if(verbose>0) printf("neg rho in u2p_solver %e %e %e %e\n",rho,uint,gamma2,W);//getchar(); - return -105; - } - + + if(rho<0.) + { + if(verbose>0) printf("neg rho in u2p_solver %e %e %e %e\n",rho,uint,gamma2,W); + //getchar(); + return -105; + } + + //entropy based on Etype //pure entropy evolution - updated only in the end of RK2 pp[ENTR]=entr; - + #ifdef MAGNFIELD //magnetic conserved=primitives pp[B1]=uu[B1] * gdetu_inv; pp[B2]=uu[B2] * gdetu_inv; pp[B3]=uu[B3] * gdetu_inv; #endif - + #ifdef EVOLVEELECTRONS conv_vels(utcon,utcon,VELPRIM,VEL4,gg,GG); ldouble Se=uu[ENTRE] * gdetu_inv / utcon[0]; pp[ENTRE]=Se; ldouble Si=uu[ENTRI] * gdetu_inv / utcon[0]; pp[ENTRI]=Si; - + #ifdef RELELECTRONS int ib; for(ib=0;ib maxerrinv) maxerrinv=errinv; - //internal energy - if(Etype==U2P_HOT) - { - iv=UU; - errinv = fabs((uunew[iv]-uu[iv])/uu[iv]); - if(errinv > maxerrinv) maxerrinv=errinv; - } - if(Etype==U2P_ENTROPY) - { - iv=ENTR; - errinv = fabs((uunew[iv]-uu[iv])/uu[iv]); - if(errinv > maxerrinv) maxerrinv=errinv; - } - - double inverr=1.e-2; - - if(Etype==U2P_ENTROPY) inverr=0.999; - - if(Etype==U2P_HOT) inverr=0.1; - - if(maxerrinv>inverr)// && verbose>0) - { - - if(Etype==U2P_ENTROPY && 0) { - printf("verify u2p (%d) failed: %e || ",Etype,maxerrinv); - printf("%e %e | %e %e | %e %e\n",uunew[RHO],uu[RHO],uunew[ENTR],uu[ENTR],uunew[VX],uu[VX]); - print_conserved(uu); - print_conserved(uunew); - print_primitives(pp); - getch(); - } - if(Etype==U2P_HOT && 0) { - printf("verify u2p (%d) failed: %e || ",Etype,maxerrinv); - printf("%e %e | %e %e | %e %e\n",uunew[RHO],uu[RHO],uunew[UU],uu[UU],uunew[VX],uu[VX]); - print_conserved(uu); - print_conserved(uunew); - print_primitives(pp); - getch(); - } - return -200; - } - } - + if(verbose>0) printf("u2p_solver returns 0\n"); return 0; //ok - } //********************************************************************** -//old Newton-Raphson solver -//iterates W, not Wp +//Newton-Raphson solver +//upgraded - uses Wp instead of W //Etype == 0 -> hot inversion (uses D,Ttt,Tti) //Etype == 1 -> entropy inversion (uses D,S,Tti) //********************************************************************** int -u2p_solver_W(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose) +u2p_solver_Wp(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose) { int i,j,k; - ldouble rho,uint,w,W,alpha,D,Sc; + ldouble rho,uint,w,W,Wp,Wpprev,alpha,D,Sc,alphasq,betasqoalphasq; ldouble ucon[4],ucov[4],utcon[4],utcov[4],ncov[4],ncon[4]; - ldouble Qcon[4],Qcov[4],Qconp[4],Qcovp[4],jmunu[4][4],Qtcon[4],Qtcov[4],Qt2,Qn; + ldouble Qcon[4],Qcov[4],Qconp[4],Qcovp[4],jmunu[4][4],Qtcon[4],Qtcov[4],Qt2,Qn,Qdotnp; ldouble QdotB,QdotBsq,Bcon[4],Bcov[4],Bsq; - - /*******************************************************/ + + /****************************/ //prepare geometry struct geometry *geom - = (struct geometry *) ggg; - + = (struct geometry *) ggg; + ldouble pgamma=GAMMA; -#ifdef CONSISTENTGAMMA + #ifdef CONSISTENTGAMMA pgamma=pick_gammagas(geom->ix,geom->iy,geom->iz); -#endif + #endif ldouble pgammam1=pgamma-1.; - + ldouble (*gg)[5], (*GG)[5], gdet, gdetu, gdetu_inv; - gg=geom->gg; GG=geom->GG; - gdet=geom->gdet; gdetu=gdet; + gg=geom->gg; GG=geom->GG; + gdet=geom->gdet;gdetu=gdet; #if (GDETIN==0) //gdet out of derivatives gdetu=1.; #endif gdetu_inv = 1. / gdetu; - + /****************************/ //equations choice int (*f_u2p)(ldouble,ldouble*,ldouble*,ldouble*,ldouble*,ldouble); - if(Etype==U2P_HOT) - f_u2p=&f_u2p_hot; - if(Etype==U2P_ENTROPY) - f_u2p=&f_u2p_entropy; + if(Etype==U2P_HOT) + f_u2p=&f_u2p_hot; + if(Etype==U2P_ENTROPY) + f_u2p=&f_u2p_entropy; /****************************/ - + if(verbose>1) { printf("********************\n"); print_conserved(uu); print_primitives(pp); } - + //conserved quantities etc //alpha - //alpha=sqrt(-1./GG[0][0]); - alpha = geom->alpha; - + alpha=geom->alpha; + alphasq=alpha*alpha; + //D D = uu[RHO] * gdetu_inv * alpha; // gdetu rho ut - - //conserved entropy "S u^t" + + //conserved entropy "gdet S u^t" Sc = uu[ENTR] * gdetu_inv * alpha; - - //Q_mu=alpha T^t_mu - Qcov[0] = (uu[UU] * gdetu_inv - uu[RHO] * gdetu_inv) * alpha; - Qcov[1] = uu[VX] * gdetu_inv * alpha; - Qcov[2] = uu[VY] * gdetu_inv * alpha; - Qcov[3] = uu[VZ] * gdetu_inv * alpha; - - //Qp_mu=alpha T^t_mu - Qcovp[0] = uu[UU] * gdetu_inv *alpha; - Qcovp[1] = uu[VX] * gdetu_inv *alpha; - Qcovp[2] = uu[VY] * gdetu_inv *alpha; - Qcovp[3] = uu[VZ] * gdetu_inv *alpha; - + + //Qp_mu=alpha T^t_mu + Qcovp[0]=uu[UU] * gdetu_inv * alpha; + Qcovp[1]=uu[VX] * gdetu_inv * alpha; + Qcovp[2]=uu[VY] * gdetu_inv * alpha; + Qcovp[3]=uu[VZ] * gdetu_inv * alpha; + //Qp^mu indices_12(Qcovp,Qconp,GG); - + + //Q_mu=alpha (T^t_mu - rho u^t delta(t,mu)) - avoid this one + Qcov[0]=(uu[UU] * gdetu_inv - uu[RHO] * gdetu_inv ) * alpha; + Qcov[1]=uu[VX] * gdetu_inv * alpha; + Qcov[2]=uu[VY] * gdetu_inv * alpha; + Qcov[3]=uu[VZ] * gdetu_inv * alpha; + //Q^mu indices_12(Qcov,Qcon,GG); - + #ifdef MAGNFIELD //curly B^mu Bcon[0]=0.; Bcon[1]=uu[B1] * gdetu_inv * alpha; Bcon[2]=uu[B2] * gdetu_inv * alpha; Bcon[3]=uu[B3] * gdetu_inv * alpha; - + //B_mu indices_21(Bcon,Bcov,gg); - + Bsq = dot(Bcon,Bcov); + QdotB = dot(Qcov,Bcon); - QdotBsq = QdotB*QdotB; + QdotBsq = QdotB*QdotB; + #else Bsq=QdotB=QdotBsq=0.; Bcon[0]=Bcon[1]=Bcon[2]=Bcon[3]=0.; -#endif // MAGNFIELD +#endif - //normal observer velocity //n_mu = (-alpha, 0, 0, 0) ncov[0]=-alpha; ncov[1]=ncov[2]=ncov[3]=0.; //n^mu indices_12(ncov,ncon,GG); - + //Q_mu n^mu = Q^mu n_mu = -alpha*Q^t - Qn = Qcon[0] * ncov[0]; + Qn=Qcon[0] * ncov[0]; + + //\beta^i \beta_i / \alpha^2 = g^{ti} g_{ti} + betasqoalphasq=gg[0][1]*GG[0][1] + gg[0][2]*GG[0][2] + gg[0][3]*GG[0][3]; - //j^mu_nu = delta^mu_nu +n^mu n_nu + //Qdotnp=-E'=-E+D + ldouble Dfactor = (-geom->gttpert + alphasq*betasqoalphasq)/(alphasq+alpha); + Qdotnp = Qconp[0]*ncov[0] + D*(Dfactor) ; // -Qdotn-W = -Qdotnp-Wp + + //j^mu_nu=delta^mu_nu +n^mu n_nu for(i=0;i<4;i++) for(j=0;j<4;j++) jmunu[i][j] = delta(i,j) + ncon[i]*ncov[j]; - + //Qtilda^nu = j^nu_mu Q^mu for(i=0;i<4;i++) { @@ -1797,22 +1435,15 @@ u2p_solver_W(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose) for(j=0;j<4;j++) Qtcon[i]+=jmunu[i][j]*Qcon[j]; } - + //Qtilda_nu indices_21(Qtcon,Qtcov,gg); - + //Qt2=Qtilda^mu Qtilda_mu Qt2=dot(Qtcon,Qtcov); FTYPE Qtsq = Qt2; - - //\beta^i \beta_i / \alpha^2 = g^{ti} g_{ti} - ldouble betasqoalphasq=gg[0][1]*GG[0][1] + gg[0][2]*GG[0][2] + gg[0][3]*GG[0][3]; - ldouble alphasq=alpha*alpha; - //Qdotnp=-E'=-E+D - ldouble Dfactor = (-geom->gttpert + alphasq*betasqoalphasq)/(alphasq+alpha); - ldouble Qdotnp = Qconp[0]*ncov[0] + D*(Dfactor) ; // -Qdotn - W = -Qdotnp-Wp - - //initial guess for W = w gamma**2 based on current primitives + + //initial guess for Wp = w gamma**2 based on primitives rho=pp[RHO]; uint=pp[UU]; utcon[0]=0.; @@ -1824,175 +1455,203 @@ u2p_solver_W(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose) { conv_vels(utcon,utcon,VELPRIM,VELR,gg,GG); } - + ldouble qsq=0.; for(i=1;i<4;i++) for(j=1;j<4;j++) qsq+=utcon[i]*utcon[j]*gg[i][j]; ldouble gamma2=1.+qsq; ldouble gamma=sqrt(gamma2); - - //W - W=(rho+pgamma*uint)*gamma2; - - if(verbose>1) printf("initial W:%e\n",W); - - // test if does not provide reasonable gamma2 - // Make sure that W is large enough so that v^2 < 1 : + + //Wp + Wp=(pgamma*uint)*gamma2; + + ldouble Wpinit, Winit; + if(verbose>1) printf("initial Wp:%e\n",Wp); + Wpinit=Wp; + Winit=Wpinit+D; + + //test if does not provide reasonable gamma2 + // Make sure that W is large enough so that v^2 < 1 and w-rho > 0 : int i_increase = 0; ldouble f0,f1,dfdW,err; - ldouble CONV=U2PCONV; - ldouble EPS=1.e-4; - ldouble Wprev=W; + ldouble CONV=U2PCONV; ldouble cons[7]={Qn,Qt2,D,QdotBsq,Bsq,Sc,Qdotnp}; + int iter=0,fu2pret; do - { - f0=dfdW=0.; - - //if(Etype!=U2P_HOT) //entropy-like solvers require this additional check - //now invoked for all solvers: - (*f_u2p)(W-D,cons,&f0,&dfdW,&err,pgamma); - - if( ((( W*W*W * ( W + 2.*Bsq ) - - QdotBsq*(2.*W + Bsq) ) <= W*W*(Qtsq-Bsq*Bsq)) - || !isfinite(f0) || !isfinite(f0) - || !isfinite(dfdW) || !isfinite(dfdW)) - && (i_increase < 50)) { - if(verbose>0) printf("init W : %e -> %e (%e %e)\n",W,100.*W,f0,dfdW); - W *= 10.; - i_increase++; - continue; + W=Wp+D; + f0=dfdW=0.; + + FTYPE Wsq,Xsq,X; + X = Bsq + W; + Xsq = X*X; + Wsq = W*W; + + ldouble v2=( Wsq * Qtsq + QdotBsq * (Bsq + 2.*W)) / (Wsq*Xsq); + ldouble gamma2 = 1./(1.-v2); + ldouble gamma = sqrt(gamma2); + ldouble rho0 = D/gamma; + ldouble wmrho0 = Wp/gamma2 - D*v2/(1.+gamma); + + (*f_u2p)(Wp,cons,&f0,&dfdW,&err,pgamma); + + if((gamma2<0. || Wp<0. || wmrho0<0.|| !isfinite(f0) || !isfinite(dfdW)) && (i_increase < 50)) + { + if(verbose>0) printf("init Wp : %e - %e %e %e %e\n",Wp,v2,wmrho0,f0,dfdW); + Wp *= 2.; + i_increase++; + continue; + } + else + break; } - else - break; - } while(1); - + if(i_increase>=50) - { - return -150; - printf("failed to find initial W for Etype: %d\n",Etype); - printf("at %d %d\n",geom->ix+TOI,geom->iy+TOJ); - print_NVvector(uu); - print_NVvector(pp); - //getchar(); - } - - //1d Newton solver - int iter=0, fu2pret; + { + if(verbose>0) + {printf("failed to find initial W for Etype: %d\n",Etype); + printf("at %d %d\n",geom->ix+TOI,geom->iy+TOJ);} + return -150; + + print_NVvector(uu); + print_NVvector(pp); + //getchar(); + } + + //1d Newton solver do - { - Wprev=W; - iter++; + { + Wpprev=Wp; + iter++; + + fu2pret=(*f_u2p)(Wp,cons,&f0,&dfdW,&err,pgamma); + + if(verbose>1) printf("%d %e %e %e %e\n",iter,Wp,f0,dfdW,err); + + //convergence test + //if(err1) printf("sub (%d) :%d %e %e %e %e %e %e\n",idump,iter,Wpnew,f0tmp,dfdWtmp,v2,gamma2,wmrho0); + + //if((gamma2<0. || Wpnew<0. || wmrho0<0. || !isfinite(f0tmp) || !isfinite(dfdWtmp)) && (idump=itmaxdamp) + { + if(verbose>0) printf("damped unsuccessfuly\n"); + return -101; + } - fu2pret=(*f_u2p)(W-D,cons,&f0,&dfdW,&err,pgamma); - - if(verbose>1) printf("%d %e %e %e %e\n",iter,W,f0,dfdW,err); - - //convergence test - if(err1) printf("sub (%d) :%d %e %e %e %e\n",idump,iter,Wnew,f0tmp,dfdWtmp,errtmp); - if( ((( Wnew*Wnew*Wnew * ( Wnew + 2.*Bsq ) - - QdotBsq*(2.*Wnew + Bsq) ) <= Wnew*Wnew*(Qtsq-Bsq*Bsq)) - || !isfinite(f0tmp) || !isfinite(f0tmp) - || !isfinite(dfdWtmp) || !isfinite(dfdWtmp)) - && (idump<100)) - { - idump++; - dumpfac/=2.; - Wnew=W-dumpfac*f0/dfdW; - continue; - } - else - break; + if(fabs(W)>BIG) + { + if(verbose>1) printf("W has gone out of bounds at %d,%d,%d\n",geom->ix+TOI,geom->iy+TOJ,geom->iz); + return -103; + } + + Wp=Wpnew; + + //convergence test: + if(err=100) + while(iter<50); + + if(iter>=50) { - if(verbose>0) printf("damped unsuccessfuly\n"); - return -101; + if(verbose>0) printf("iter exceeded in u2p_solver with Etype: %d\n",Etype); //getchar(); + return -102; } - - W=Wnew; - - if(fabs(W)>BIG) + + if(!isfinite(Wp) || !isfinite(Wp)) {if(verbose) printf("nan/inf W in u2p_solver with Etype: %d\n",Etype); return -103;} + + if(verbose>1) { - if(verbose>1) printf("W has gone out of bounds at %d,%d,%d\n",geom->ix+TOI,geom->iy+TOJ,geom->iz); - return -103; + fu2pret=(*f_u2p)(Wp,cons,&f0,&dfdW,&err,pgamma); + printf("end: %d %e %e %e %e\n",iter,Wp,f0,dfdW,err); } - - if(fabs((W-Wprev)/Wprev)=50) - { - if(verbose>0) printf("iter exceeded in u2p_solver with Etype: %d\n",Etype); //getchar(); - return -102; - } - - if(!isfinite(W)) - { - if(verbose) printf("nan/inf W in u2p_solver with Etype: %d\n",Etype); - return -103; - } - - if(verbose>1) - { - fu2pret=(*f_u2p)(W-D,cons,&f0,&dfdW,&err,pgamma); - printf("end: %d %e %e %e %e\n",iter,W,f0,dfdW,err); - } - + //W found, let's calculate v2 and the rest - ldouble Wsq,Xsq,v2,entr; - + W=Wp+D; + ldouble Wsq,Xsq,v2,wmrho0,entr; + Wsq = W*W ; - Xsq = (Bsq + W) * (Bsq + W); + Xsq = (Bsq + W) * (Bsq + W); v2 = ( Wsq * Qtsq + QdotBsq * (Bsq + 2.*W)) / (Wsq*Xsq); - + gamma2=1./(1.-v2); gamma=sqrt(gamma2); rho=D/gamma; entr=Sc/gamma; - uint=1./pgamma*(W/gamma2-rho); + // w-\rho_0 = (u+p) = W'/\gamma^2 - D v^2/(1+\gamma) + wmrho0 = Wp/gamma2 - D*v2/(1.+gamma); + uint=1./pgamma*wmrho0; utcon[0]=0.; utcon[1]=gamma/(W+Bsq)*(Qtcon[1]+QdotB*Bcon[1]/W); utcon[2]=gamma/(W+Bsq)*(Qtcon[2]+QdotB*Bcon[2]/W); utcon[3]=gamma/(W+Bsq)*(Qtcon[3]+QdotB*Bcon[3]/W); - if(verbose>1) - printf("end2: %e %e %e %e %e %e\n",W,D,pgamma,gamma2,rho,uint); if(!isfinite(utcon[1])) - { - return -120; - } - - if(uint<0. || gamma2<0. || isnan(W) || !isfinite(W)) - { - if(verbose>0) printf("neg u in u2p_solver %e %e %e %e\n",rho,uint,gamma2,W);//getchar(); - return -104; - } - + { + //print_4vector(utcon); + return -120; + } + + + if(uint<0. || gamma2<0. ||isnan(Wp) || !isfinite(Wp)) + { + if(verbose>0) printf("neg u in u2p_solver %e %e %e %e\n",rho,uint,gamma2,W);//getchar(); + return -104; + } + //converting to VELPRIM if (VELR != VELPRIM) { @@ -2001,48 +1660,207 @@ u2p_solver_W(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose) //returning new primitives pp[RHO]=rho; - pp[UU]=uint; - pp[VX]=utcon[1]; - pp[VY]=utcon[2]; - pp[VZ]=utcon[3]; - - if(rho<0.) - { - if(verbose>0) printf("neg rho in u2p_solver %e %e %e %e\n",rho,uint,gamma2,W); - //getchar(); - return -105; - } - - //entropy based on Etype - //pure entropy evolution - updated only in the end of RK2 - pp[ENTR]=entr; - + pp[UU]=uint; + pp[VX]=utcon[1]; + pp[VY]=utcon[2]; + pp[VZ]=utcon[3]; + + if(rho<0.) + { + if(verbose>0) printf("neg rho in u2p_solver %e %e %e %e\n",rho,uint,gamma2,W);//getchar(); + return -105; + } + + //pure entropy evolution - updated only in the end of RK2 + pp[ENTR]=entr; + +#ifdef MAGNFIELD + //magnetic conserved=primitives + pp[B1]=uu[B1] * gdetu_inv; + pp[B2]=uu[B2] * gdetu_inv; + pp[B3]=uu[B3] * gdetu_inv; +#endif + +#ifdef EVOLVEELECTRONS + conv_vels(utcon,utcon,VELPRIM,VEL4,gg,GG); + ldouble Se=uu[ENTRE] * gdetu_inv / utcon[0]; + pp[ENTRE]=Se; + ldouble Si=uu[ENTRI] * gdetu_inv / utcon[0]; + pp[ENTRI]=Si; + +#ifdef RELELECTRONS + int ib; + for(ib=0;ib maxerrinv) maxerrinv=errinv; + //internal energy + if(Etype==U2P_HOT) + { + iv=UU; + errinv = fabs((uunew[iv]-uu[iv])/uu[iv]); + if(errinv > maxerrinv) maxerrinv=errinv; + } + if(Etype==U2P_ENTROPY) + { + iv=ENTR; + errinv = fabs((uunew[iv]-uu[iv])/uu[iv]); + if(errinv > maxerrinv) maxerrinv=errinv; + } + + double inverr=1.e-2; + + if(Etype==U2P_ENTROPY) inverr=0.999; + + if(Etype==U2P_HOT) inverr=0.1; + + if(maxerrinv>inverr)// && verbose>0) + { + + if(Etype==U2P_ENTROPY && 0) { + printf("verify u2p (%d) failed: %e || ",Etype,maxerrinv); + printf("%e %e | %e %e | %e %e\n",uunew[RHO],uu[RHO],uunew[ENTR],uu[ENTR],uunew[VX],uu[VX]); + print_conserved(uu); + print_conserved(uunew); + print_primitives(pp); + getch(); + } + if(Etype==U2P_HOT && 0) { + printf("verify u2p (%d) failed: %e || ",Etype,maxerrinv); + printf("%e %e | %e %e | %e %e\n",uunew[RHO],uu[RHO],uunew[UU],uu[UU],uunew[VX],uu[VX]); + print_conserved(uu); + print_conserved(uunew); + print_primitives(pp); + getch(); + } + return -200; + } + } + + if(verbose>0) printf("u2p_solver returns 0\n"); + return 0; //ok + +} + +//********************************************************************** +//non-relativistic, analytical, solver +//********************************************************************** + +int +u2p_solver_nonrel(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose) +{ + //prepare geometry + struct geometry *geom + = (struct geometry *) ggg; + + ldouble (*gg)[5], (*GG)[5], gdet, gdetu, gdetu_inv; + gg=geom->gg; GG=geom->GG; + gdet=geom->gdet;gdetu=gdet; +#if (GDETIN==0) //gdet out of derivatives + gdetu=1.; +#endif + gdetu_inv = 1. / gdetu; + + //density + ldouble rho=uu[RHO] * gdetu_inv; + pp[RHO]=rho; + + //velocities u_i + ldouble ucov[4],ucon[4],vcov[4]; + ucov[0]=-1.; + ucov[1]=uu[VX] * gdetu_inv / rho; + ucov[2]=uu[VY] * gdetu_inv / rho; + ucov[3]=uu[VZ] * gdetu_inv / rho; + + indices_12(ucov,ucon,GG); + + ucon[0]=1.; + + ldouble v2=dot3nr(ucon,ucov); + + pp[VX]=ucon[1]; + pp[VY]=ucon[2]; + pp[VZ]=ucon[3]; + + ldouble bsq=0.; + #ifdef MAGNFIELD - //magnetic conserved=primitives - pp[B1]=uu[B1] * gdetu_inv; - pp[B2]=uu[B2] * gdetu_inv; - pp[B3]=uu[B3] * gdetu_inv; + + ldouble bcon[4],bcov[4],Bcon[4]; + + Bcon[0]=0.; + Bcon[1]=uu[B1] * gdetu_inv; + Bcon[2]=uu[B2] * gdetu_inv ; + Bcon[3]=uu[B3] * gdetu_inv ; + + pp[B1]=Bcon[1]; + pp[B2]=Bcon[2]; + pp[B3]=Bcon[3]; + + calc_bcon_bcov_bsq_from_4vel(pp, ucon, ucov, geom, bcon, bcov, &bsq); + #endif + -#ifdef EVOLVEELECTRONS - conv_vels(utcon,utcon,VELPRIM,VEL4,gg,GG); - ldouble Se=uu[ENTRE] * gdetu_inv / utcon[0]; + ldouble uint; + if(Etype==U2P_HOT) + { + uint = -uu[UU] * gdetu_inv - bsq/2. - rho*v2/2.; + + if(uint %e %e %e %e %e\n",geom->ix+TOI,geom->iy+TOI,uint,uu[UU]/gdetu,bsq/2.,rho*v2/2.,rho); + //printf("%e %e\n",uint,rho); + //if(geom->ix>50) getch(); + return -1; + } + + pp[UU]=uint; + } + else if(Etype==U2P_ENTROPY) + { + ldouble S=uu[ENTR] * gdetu_inv ; + uint= calc_ufromS(S,rho,geom->ix,geom->iy,geom->iz); + pp[UU]=uint; + } + + //pure entropy evolution - updated only in the end of RK2 + pp[ENTR]= uu[ENTR] * gdetu_inv; + + #ifdef EVOLVEELECTRONS + ldouble Se=uu[ENTRE] * gdetu_inv ; pp[ENTRE]=Se; - ldouble Si=uu[ENTRI] * gdetu_inv / utcon[0]; + + ldouble Si=uu[ENTRI] * gdetu_inv ; pp[ENTRI]=Si; - -#ifdef RELELECTRONS - int ib; - for(ib=0;ib0) printf("u2p_solver returns 0\n"); - return 0; //ok -} + #endif + + #ifdef RELELECTRONS + int ib; + for(ib=0;ibgg; GG=geom->GG; - gdet=geom->gdet;gdetu=gdet; + FTYPE W3,X3,Ssq,Wsq,X,X2,Xsq; + FTYPE Qtsq = Qt2; + X = Bsq + W; + Wsq = W*W; + W3 = Wsq*W ; + X2 = X*X; + Xsq = X2; + X3 = X2*X; - #if (GDETIN==0) //gdet out of derivatives - gdetu=1.; - #endif - gdetu_inv = 1. / gdetu; + ldouble v2=( Wsq * Qtsq + QdotBsq * (Bsq + 2.*W)) / (Wsq*Xsq); + ldouble gamma2 = 1./(1.-v2); + ldouble gamma = sqrt(gamma2); + ldouble rho0 = D/gamma; + ldouble wmrho0 = Wp/gamma2 - D*v2/(1.+gamma); + ldouble u = wmrho0 / pgamma; + ldouble p = (pgamma-1)*u; - ldouble pgamma=GAMMA; - #ifdef CONSISTENTGAMMA - pgamma=pick_gammagas(geom->ix,geom->iy,geom->iz); - #endif + //original: +#if (U2P_EQS==U2P_EQS_NOBLE) + *f = Qn + W - p + 0.5*Bsq*(1.+v2) - QdotBsq/2./Wsq; + *err = fabs(*f) / (fabs(Qn) + fabs(W) + fabs(p) + fabs(0.5*Bsq*(1.+v2)) + fabs(QdotBsq/2./Wsq)); +#endif - //alpha - alpha = geom->alpha; - alphasq = alpha*alpha; + //JONS: +#if (U2P_EQS==U2P_EQS_JON) + *f = Qdotnp + Wp - p + 0.5*Bsq + (Bsq*Qtsq - QdotBsq)/X2; + *err = fabs(*f) / (fabs(Qdotnp) + fabs(Wp) + fabs(p) + fabs(0.5*Bsq) + fabs((Bsq*Qtsq - QdotBsq)/X2)); +#endif - //beta vector components - betavec[0] = 0.; - betavec[1] = alphasq * GG[0][1]; - betavec[2] = alphasq * GG[0][2]; - betavec[3] = alphasq * GG[0][3]; + // dp/dWp = dp/dW + dP/dv^2 dv^2/dW + ldouble dvsq=(-2.0/X3 * ( Qtsq + QdotBsq * (3.0*W*X + Bsq*Bsq)/W3)); + ldouble dp1 = dpdWp_calc_vsq(Wp, D, v2 ,pgamma); // vsq can be unphysical - ////////////////////////////////////////////////////////////////////////////////////////// - // Force-Free EM part - ///////////////////////////////////////////////////////////////////////////////////////// + ldouble idwmrho0dp=compute_idwmrho0dp(wmrho0,pgamma); + ldouble dwmrho0dvsq = (D*(gamma*0.5-1.0) - Wp); - //McKinney defn of B^mu = B^mu HARM * alpha - Bcon[0]=0.; - Bcon[1]=uu[B1] * alpha * gdetu_inv; - Bcon[2]=uu[B2] * alpha * gdetu_inv; - Bcon[3]=uu[B3] * alpha * gdetu_inv; + ldouble drho0dvsq = -D*gamma*0.5; // because \rho = D/\gamma + ldouble idrho0dp = compute_idrho0dp(wmrho0); - indices_21(Bcon,Bcov,gg); - Bsq = dot(Bcon,Bcov); - Bmag = sqrt(Bsq); - - //Q_mu= alpha*T^t_mu - Qcov[0] = 0.; // this is not the actual Qcov[0] for FF, since we only conserve spatial comps - //Qcov[0] = uu[UU] * alpha * gdetu_inv; - Qcov[1] = uu[VXFF] * alpha * gdetu_inv; - Qcov[2] = uu[VYFF] * alpha * gdetu_inv; - Qcov[3] = uu[VZFF] * alpha * gdetu_inv; - - //Qtilde_mu = Q_mu + n_mu (n.Q) - Qtildecov[0] = Qcov[1]*betavec[1] + Qcov[2]*betavec[2] + Qcov[3]*betavec[3]; - Qtildecov[1] = Qcov[1]; - Qtildecov[2] = Qcov[2]; - Qtildecov[3] = Qcov[3]; - indices_12(Qtildecov,Qtildecon,GG); - - // get three velocity - vcon_perp[0] = 0.; - vcon_perp[1] = Qtildecon[1]/Bsq; - vcon_perp[2] = Qtildecon[2]/Bsq; - vcon_perp[3] = Qtildecon[3]/Bsq; - - // get Lorentz factor - int ll,mm; - ldouble vsq_perp=0.; - ldouble vsqmax = 1. - 1./(GAMMAMAXFF*GAMMAMAXFF); - int hitgammaceil=0; - for(ll=1;ll<4;ll++) - { - for(mm=1;mm<4;mm++) - { - vsq_perp += gg[ll][mm]*vcon_perp[ll]*vcon_perp[mm]; - } - } + ldouble dp2 = drho0dvsq *idrho0dp + dwmrho0dvsq *idwmrho0dp; - ldouble vsq_perp_orig = vsq_perp; + ldouble dpdW = dp1 + dp2*dvsq; // dp/dW = dp/dWp - // lorentz factor limiter - if(vsq_perp>vsqmax) - { - hitgammaceil=1; - ldouble vfac = sqrt(vsqmax/vsq_perp); - vcon_perp[1] *= vfac; - vcon_perp[2] *= vfac; - vcon_perp[3] *= vfac; - vsq_perp = vsqmax; - } + //original: + #if (U2P_EQS==U2P_EQS_NOBLE) + *df=1. -dpdW + QdotBsq/(Wsq*W) + 0.5*Bsq*dvsq; + #endif - gamma2_perp = 1./(1-vsq_perp); - gamma_perp = sqrt(gamma2_perp); - - if(hitgammaceil==1) - { - //ldouble alpgam = calc_alpgam(utcon,gg,GG); - //printf("LIMITER vsq orig: %e | new: %e | gammma %e \n",vsq_orig,vsq,gamma); getch(); - } - - if(!isfinite(gamma_perp)) - { - if(verbose>0) printf("gamma_perp not finite in u2p_solver_ff %e \n",Bsq); - //getchar(); - return -104; - } + //JONS: + #if (U2P_EQS==U2P_EQS_JON) + *df=1. -dpdW + (Bsq*Qtsq - QdotBsq)/X3*(-2.0); + #endif + return 0; +} - //************************************** - // Return new primitives - pp[VXFF]=vcon_perp[1]*gamma_perp; - pp[VYFF]=vcon_perp[2]*gamma_perp; - pp[VZFF]=vcon_perp[3]*gamma_perp; +//******************************************** +//Harm u2p_entropy +//******************************************** - pp[B1]=uu[B1] * gdetu_inv; // TODO or Bcon[1]/alpha; - pp[B2]=uu[B2] * gdetu_inv; - pp[B3]=uu[B3] * gdetu_inv; +// p(rho0, w-rho0 = u+p) +static FTYPE +pressure_wmrho0_idealgas(FTYPE rho0, FTYPE wmrho0,FTYPE gamma) +{ + ldouble igammar = (gamma-1.)/gamma; + return(igammar*wmrho0) ; +} + +// local aux function +static FTYPE +compute_inside_entropy_wmrho0_idealgas(FTYPE rho0, FTYPE wmrho0,FTYPE gamma) +{ + FTYPE pressure,indexn,insideentropy; - ////////////////////////////////////////////////////////////////////////////////////////// - // Hydro part -- solve for parallel velocity - ///////////////////////////////////////////////////////////////////////////////////////// - - ldouble gammam2_perp=1./gamma2_perp; - ldouble afac=(pgamma-1.)/pgamma; - - // MHD conserveds - ldouble D = uu[RHO] * alpha * gdetu_inv; // gdet rho ut, so D = gamma * rho - ldouble Sc = uu[ENTR] * alpha * gdetu_inv; // gdet S ut, so Sc = gamma * S - ldouble W = -1.; - ldouble Y=0; - ldouble Z=0; - -#ifdef FORCEFREE_SOLVE_PARALLEL - Qcov[0] = uu[UU] * alpha * gdetu_inv; - Qcov[1] = uu[VX] * alpha * gdetu_inv; - Qcov[2] = uu[VY] * alpha * gdetu_inv; - Qcov[3] = uu[VZ] * alpha * gdetu_inv; - - indices_12(Qcov,Qcon,GG); - - ldouble QdotB = Bcon[1]*Qcov[1] + Bcon[2]*Qcov[2] + Bcon[3]*Qcov[3]; // TODO verify - ldouble QdotEta = -alpha * Qcon[0]; - - Y = QdotB / Bmag; - Z = QdotEta + 0.5*Bsq*vsq_perp; - - ldouble cons[5]; - cons[0] = D; - cons[1] = Y; - cons[2] = Z; - cons[3] = gammam2_perp; - cons[4] = afac; - - // initial guess is based on current primitives TODO - // note force free requires VELPRIM = VELR - ldouble utcon_prev[4]; - utcon_prev[0]=0.; - utcon_prev[1]=pp[VX]; - utcon_prev[2]=pp[VY]; - utcon_prev[3]=pp[VZ]; - - ldouble qsq=0.; - for(i=1;i<4;i++) - for(j=1;j<4;j++) - qsq+=utcon_prev[i]*utcon_prev[j]*gg[i][j]; - ldouble gamma2_prev=1.+qsq; + pressure=pressure_wmrho0_idealgas(rho0,wmrho0,gamma); + indexn=1.0/(gamma-1.); + insideentropy=pow(pressure,indexn)/pow(rho0,indexn+1.0); - Wguess = (pp[RHO] + pgamma*pp[UU])*gamma2_prev; + return(insideentropy); +} - // solve for W - W = u2p_solver_ff_parallel(Wguess,cons,1); -#endif // FORCEFREE_SOLVE_PARALLEL - - // acceptable solution on parallel momentum solver - // TODO other conditions? - if(W>0) - { - // parallel three-velocity - ldouble vpar = Y / W; - - vcon_par[0] = 0.; - vcon_par[1] = vpar*(Bcon[1]/Bmag); - vcon_par[2] = vpar*(Bcon[2]/Bmag); - vcon_par[3] = vpar*(Bcon[3]/Bmag); - - // get parallel Lorentz factor - ldouble vsq_par=0.; - vsqmax = 1. - 1./(GAMMAMAXFF*GAMMAMAXFF); - hitgammaceil=0; - for(ll=1;ll<4;ll++) - { - for(mm=1;mm<4;mm++) - { - vsq_par += gg[ll][mm]*vcon_par[ll]*vcon_par[mm]; - } - } - // limiter - if(vsq_par>vsqmax) - { - hitgammaceil=1; - ldouble vfac = sqrt(vsqmax/vsq_par); - vcon_par[1] *= vfac; - vcon_par[2] *= vfac; - vcon_par[3] *= vfac; - vsq_par = vsqmax; - } - - gamma2_par = 1./(1-vsq_par); - gamma_par = sqrt(gamma2_par); - - gamma2 = 1./(1 - vsq_perp - vsq_par); - gamma = sqrt(gamma2); - - // total velocity (VELR) - pp[VX]=(vcon_perp[1] + vcon_par[1])*gamma; - pp[VY]=(vcon_perp[2] + vcon_par[2])*gamma; - pp[VZ]=(vcon_perp[3] + vcon_par[3])*gamma; - - // density, entropy, energy density - rho = D/gamma; - entr = Sc/gamma; - uint = (W/gamma2 - rho)/pgamma; - - } +// specific entropy as function of rho0 and internal energy (u) +// Ss(rho0,\chi=u+p) +// specific entropy = \ln( p^n/\rho^{n+1} ) +static FTYPE +compute_specificentropy_wmrho0_idealgas(FTYPE rho0, FTYPE wmrho0,FTYPE gamma) +{ + FTYPE insideentropy,specificentropy; - // unacceptable solution, do adiabatic with NO parallel velocity - else //(W<0) - { - if(verbose>1) printf("reverting to adiabatic in FF+MHD solver!\n"); + insideentropy=compute_inside_entropy_wmrho0_idealgas(rho0, wmrho0,gamma); + specificentropy=log(insideentropy); - // total velocity (VELR) has no parallel component - pp[VX] = pp[VXFF]; - pp[VY] = pp[VYFF]; - pp[VZ] = pp[VZFF]; + return(specificentropy); - rho = D/gamma_perp; - entr = Sc/gamma_perp; - uint = calc_ufromS(entr,rho,geom->ix,geom->iy,geom->iz); - } - - if(rho<0. || !isfinite(rho)) - { - if(verbose>0) printf("neg rho in u2p_solver_ff %e %e %e \n",rho,D,gamma); //getchar(); - return -105; - } - - if(uint<0. || !isfinite(uint)) - { - if(verbose>0) printf("neg uint in u2p_solver_ff %e %e %e\n",uint,entr,rho); //getchar(); - return -107; - } - - //*************************************** - //final new primitives - pp[RHO] = rho; - pp[UU] = uint; - pp[ENTR] = entr; - - //if(verbose>0) printf("u2p_solver_ff returns 0\n"); - return 0; //ok -#endif } -//********************************************************************* -//parallel momentum part of the force-free solver -//********************************************************************* -// TODO solve for W? Wp? vparallel^2? -// TODO exact solution? -// TODO entropy version? -// TODO what to do about variable adiabatic index? -static ldouble -u2p_solver_ff_parallel(ldouble Wguess,ldouble* cons, int verbose) +// used for utoprim_jon when doing entropy evolution +// dSspecific/d\chi +static FTYPE +compute_dspecificSdwmrho0_wmrho0_idealgas(FTYPE rho0, FTYPE wmrho0,FTYPE gamma) { - ldouble CONV=U2PCONV; - ldouble EPS=1.e-4; - - ldouble W=Wguess; - - if(verbose>1) printf("initial W:%e\n",W); - - // constants - ldouble D = cons[0]; - ldouble Y = cons[1]; - ldouble Z = cons[2]; - ldouble gammam2_perp = cons[3]; - ldouble afac = cons[4]; - - ldouble Ysq=Y*Y; - - ldouble f0,f1,dfdW,err,vpar2guess; - ldouble Wprev=W; - int i_increase = 0; - // Make sure that W is large enough so that v^2 < 1 : - // Test if initial guess is out of bounds and damp if so - do - { - f0=dfdW=0.; - - (*f_u2p_parallel)(W,cons,&f0,&dfdW,&err); + FTYPE dSdchi; - vpar2guess = Ysq/(W*W); - if(((gammam2_perp - vpar2guess) < 0 - || !isfinite(f0) || !isfinite(dfdW) || !isfinite(dfdW)) - && (i_increase < 50)) - { - if(verbose>0) printf("init W : %e -> %e (%e %e)\n",W,100.*W,f0,dfdW); - W *= 10.; - i_increase++; - continue; - } - else - break; - } - while(1); - - if(i_increase>=50) - { - printf("failed to find initial W in parallel solver"); - //printf("at %d %d\n",geom->ix+TOI,geom->iy+TOJ); - //print_NVvector(uu); - //print_NVvector(pp); - getchar(); - return -150; - } - - //1d Newton solver - int iter=0, fu2pret; - do - { - Wprev=W; - iter++; + dSdchi = 1.0/((gamma-1.)*wmrho0); + // Again, GAMMA->1 means dSdchi->\infty unless \chi->0 or rho0->0 - fu2pret=(*f_u2p_parallel)(W,cons,&f0,&dfdW,&err); - - if(verbose>1) printf("parallel solver: %d %e %e %e %e\n",iter,W,f0,dfdW,err); - - //convergence test - if(err=100) - { - if(verbose>0) printf("damped unsuccessfuly in parallel solver\n"); - getchar(); - return -101; - } - - W=Wnew; - - if(fabs(W)>BIG) - { - if(verbose>1) printf("W has gone out of bounds in parallel solver\n"); - getchar(); - return -103; - } +} - // TODO right error? - if(fabs((W-Wprev)/Wprev)=50) - { - if(verbose>0) printf("iter exceeded in parallel u2p_solver \n"); - getchar(); - return -102; - } +// dSspecific/drho0 +static FTYPE +compute_dspecificSdrho_wmrho0_idealgas(FTYPE rho0, FTYPE wmrho0, FTYPE gamma) +{ + FTYPE dSdrho; - if(!isfinite(W)) - { - if(verbose) printf("nan/inf W in parallel u2p_solver: %d\n"); - getchar(); - return -103; - } - - vpar2guess = Ysq/(W*W); - if((gammam2_perp - vpar2guess) < 0) - { - if(verbose>0) printf("final W out of bounds in parallel solver\n"); - getchar(); - return -109; - } + dSdrho=gamma/((1.0-gamma)*rho0); - return W; + return(dSdrho); } -// TODO variable adiabatic index? ? static int -f_u2p_parallel(ldouble W, ldouble* cons,ldouble *f,ldouble *df,ldouble *err) +f_u2p_entropy(ldouble Wp, ldouble* cons, ldouble *f, ldouble *df, ldouble *err,ldouble pgamma) { + ldouble Qn=cons[0]; + ldouble Qt2=cons[1]; + ldouble D=cons[2]; + ldouble QdotBsq=cons[3]; + ldouble Bsq=cons[4]; + ldouble Sc=cons[5]; + + ldouble W=Wp+D; - // constants - ldouble D = cons[0]; - ldouble Y = cons[1]; - ldouble Z = cons[2]; - ldouble gammam2_perp = cons[3]; - ldouble afac = cons[4]; + FTYPE W3,X3,Ssq,Wsq,X,X2,Xsq; + FTYPE Qtsq = Qt2; + X = Bsq + W; + Wsq = W*W; + W3 = Wsq*W ; + X2 = X*X; + Xsq = X2; + X3 = X2*X; - ldouble Ysq = Y*Y; - ldouble Wsq = W*W; + ldouble v2=( Wsq * Qtsq + QdotBsq * (Bsq + 2.*W)) / (Wsq*Xsq); + ldouble gamma2 = 1./(1.-v2); + ldouble gamma = sqrt(gamma2); + ldouble rho0 = D/gamma; + ldouble wmrho0 = Wp/gamma2 - D*v2/(1.+gamma); + ldouble u = wmrho0 / pgamma; + ldouble p = (pgamma-1)*u; - - // root function - *f = afac*W*(gammam2_perp - Ysq/Wsq) - afac*D*sqrt(gammam2_perp - Ysq/Wsq) - W - Z; + ldouble Ssofchi=compute_specificentropy_wmrho0_idealgas(rho0,wmrho0,pgamma); - // error - *err = fabs(*f) / (fabs(afac*W*(gammam2_perp - Ysq/Wsq)) - fabs(afac*D*sqrt(gammam2_perp - Ysq/Wsq)) - fabs(W) - fabs(Z)); + *f= -Sc + D*Ssofchi; + *err = fabs(*f) / (fabs(Sc) + fabs(D*Ssofchi)); - // derivitive - ldouble dterm1dW = afac*(gammam2_perp + Ysq/Wsq); - ldouble dterm2dW = -afac*D*Ysq / (W*W*W*sqrt(gammam2_perp - Ysq/Wsq)); - ldouble dterm3dW = -1.; - ldouble dterm4dW = 0.; + FTYPE dSsdW,dSsdvsq,dSsdWp,dScprimedWp,dSsdrho,dSsdchi; + FTYPE dvsq,dwmrho0dW,drho0dW; + FTYPE dwmrho0dvsq,drho0dvsq; - *df = dterm1dW + dterm2dW + dterm3dW + dterm4dW; + dSsdrho=compute_dspecificSdrho_wmrho0_idealgas(rho0,wmrho0,pgamma); + dSsdchi=compute_dspecificSdwmrho0_wmrho0_idealgas(rho0,wmrho0,pgamma); - return 0; -} + dwmrho0dW = 1.0/gamma; // holding utsq fixed + drho0dW = 0.0; // because \rho=D/\gamma and holding utsq fixed + dwmrho0dvsq = (D*(gamma*0.5-1.0) - Wp); // holding Wp fixed + drho0dvsq = -D*gamma*0.5; // because \rho=D/\gamma and holding Wp fixed + + dvsq=(-2.0/X3 * ( Qtsq + QdotBsq * (3.0*W*X + Bsq*Bsq)/W3)); + + dSsdW = drho0dW *dSsdrho + dwmrho0dW *dSsdchi; // dSs/dW' holding utsq fixed + dSsdvsq = drho0dvsq *dSsdrho + dwmrho0dvsq *dSsdchi; + dSsdWp = dSsdW + dSsdvsq*dvsq; // dSs/dW = dSs/dWp [total derivative] + + dScprimedWp = D*dSsdWp; + *df = dScprimedWp; + + return 0; + +} //********************************************************************** //count the number of entropy inversions @@ -2645,37 +2214,148 @@ update_entropy() int test_inversion() { - ldouble pp[NV],pp2[NV],uu[NV],ucon[4]={0.,0.,0.,0.}; + ldouble pp[NV],pp0[NV],uu[NV],vcon[4]={0.,0.,0.,0.},Bcon[4]={0.,0.,0.,0.}; + ldouble Bcov[4]; struct geometry geom,geomBL; int iv; - fill_geometry(NX-2,0,0,&geom); - fill_geometry_arb(NX-2,0,0,&geomBL,BLCOORDS); - ucon[1]=1.e-12; - conv_vels(ucon,ucon,VEL4,VEL4,geomBL.gg,geomBL.GG); - trans2_coco(geomBL.xxvec,ucon,ucon,BLCOORDS,MYCOORDS); - conv_vels(ucon,ucon,VEL4,VELPRIM,geom.gg,geom.GG); - pp[RHO]=1.; - pp[UU]=1.; - pp[VX]=ucon[1]; - pp[VY]=pp[VZ]=0.; + int ix=10, iy=8, iz=0; + fill_geometry(ix,iy,iz,&geom); + ldouble alpha=geom.alpha; + ldouble (*gg)[5]; + gg=geom.gg; + //fill_geometry_arb(ix,iy,iz,&geomBL,BLCOORDS); + + ldouble rho=1.3; + ldouble uint=0.424; + + vcon[0]=0.; + vcon[1]=1.e-2; + vcon[2]=5.e-3; + vcon[3]=3.e-5; + + Bcon[0]=0.; + Bcon[1]=10.*alpha; + Bcon[2]=1.e-2*alpha; + Bcon[3]=1.e-3*alpha; + + indices_21(Bcon,Bcov,gg); + + ldouble Bsq=dot(Bcon,Bcov); + ldouble vpar=dot(vcon,Bcov)/sqrt(Bsq); + + ldouble vcon_perp[4],vcon_par[4]; + vcon_par[0] = 0.; + vcon_par[1] = vpar*Bcon[1]/sqrt(Bsq); + vcon_par[2] = vpar*Bcon[2]/sqrt(Bsq); + vcon_par[3] = vpar*Bcon[3]/sqrt(Bsq); + vcon_perp[0] = 0.; + vcon_perp[1] = vcon[1]-vcon_par[1]; + vcon_perp[2] = vcon[2]-vcon_par[2]; + vcon_perp[3] = vcon[3]-vcon_par[3]; + + printf("perptest! %e %e %e\n", dot(Bcov,vcon),dot(Bcov,vcon_par),dot(Bcov,vcon_perp)); + + ldouble vsq=0.,vsq_perp=0.,vsq_par=0.; + int ll,mm; + for(ll=1;ll<4;ll++) + { + for(mm=1;mm<4;mm++) + { + vsq += gg[ll][mm]*vcon[ll]*vcon[mm]; + vsq_perp += gg[ll][mm]*vcon_perp[ll]*vcon_perp[mm]; + vsq_par += gg[ll][mm]*vcon_par[ll]*vcon_par[mm]; + } + } + + ldouble gamma = sqrt(1./(1-vsq)); + ldouble gamma_perp = sqrt(1./(1-vsq_perp)); + ldouble gamma_par = sqrt(1./(1-vsq_par)); + printf("gamma %e , gamma_perp %e, gamma_par %e, check %e\n",gamma,gamma_perp,gamma_par,1./(gamma_perp*gamma_perp)+1./(gamma_par*gamma_par)-1./(gamma*gamma)); + + ldouble afac=(GAMMA-1.)/GAMMA; + ldouble p = (GAMMA-1.)*uint; + ldouble D = rho*gamma; + ldouble W = (rho+GAMMA*uint)*gamma*gamma; + ldouble Y = W*vpar; + ldouble QdotEta = -0.5*Bsq*(1+vsq) - W + p + 0.5*vsq_par*Bsq; + ldouble Z = QdotEta + 0.5*Bsq*(1+vsq_perp); + + ldouble Sc = gamma*calc_Sfromu(rho,uint,ix,iy,iz); + + printf("sigma %e\n",Bsq/(gamma_perp*gamma_perp)/rho); + printf("Bsq %e\n",Bsq); + printf("W %e\n",W); + printf("D %e Y %e Z %e\n",D,Y,Z); + printf("gammam2perp %e\n",1./gamma_perp/gamma_perp); + printf("Sc %e \n",Sc); + printf("check1 %e\n",Z-(p-W)); + printf("check2 %e\n",(afac/(gamma*gamma)-1)*W - afac*D/gamma - Z); + + // put primitives in + pp[RHO]=rho; + pp[UU]=uint; + pp[ENTR] = calc_Sfromu(rho,uint,ix,iy,iz); + + pp[VX]=gamma*vcon[1]; + pp[VY]=gamma*vcon[2]; + pp[VZ]=gamma*vcon[3]; + + pp[B1] = Bcon[1]/alpha; + pp[B2] = Bcon[2]/alpha; + pp[B3] = Bcon[3]/alpha; + + #ifdef FORCEFREE + pp[UUFF] = uint; + pp[VXFF] = gamma_perp*vcon_perp[1]; + pp[VYFF] = gamma_perp*vcon_perp[2]; + pp[VZFF] = gamma_perp*vcon_perp[3]; + #endif + + printf("input\n"); print_primitives(pp); p2u(pp,uu,&geom); - pp[VX]*=100000.*M_PI; - pp[RHO]*=0.001245325124; + PLOOP(iv) pp0[iv]=pp[iv]; + + pp[RHO]*=3.1245325124; pp[UU]*=23.124124214421124; + //pp[ENTR] = calc_Sfromu(pp[RHO],pp[UU],ix,iy,iz); - PLOOP(iv) pp2[iv]=pp[iv]; + pp[VX]*=12.3; + pp[VY]*=21.1; + pp[VZ]*=.66; - - print_conserved(uu); - printf("gdet = %e\n",geom.gdet); - u2p_solver(uu,pp,&geom,U2P_HOT,0); + #ifdef FORCEFREE + pp[UUFF]=pp[UU]; + pp[VXFF]*=22.3; + pp[VYFF]*=12.1; + pp[VZFF]*=5.6; + #endif + + //print_conserved(uu); + u2p_solver_ff(uu,pp,&geom,U2P_ENTROPY,2); + //u2p_solver(uu,pp,&geom,U2P_HOT,0); + + printf("solved\n"); print_primitives(pp); + ldouble pperr[NV]; + PLOOP(iv) pperr[iv]=(pp[iv]-pp0[iv])/pp0[iv]; + + printf("error\n"); + print_primitives(pperr); + + printf("projection test....\n"); + PLOOP(iv) pp[iv]=pp0[iv]; + update_ffprims_cell(pp, uu, &geom); + print_primitives(pp); + PLOOP(iv) pperr[iv]=(pp[iv]-pp0[iv])/pp0[iv]; + printf("projection error\n"); + print_primitives(pperr); + return 0; } diff --git a/u2prad.c b/u2prad.c index 6ccef6d..5819ce6 100644 --- a/u2prad.c +++ b/u2prad.c @@ -532,6 +532,7 @@ check_floors_rad(ldouble *pp, int whichvel,void *ggg) int verbose=0; int ret=0; +#if !defined(SKIPALLFLOORS) struct geometry *geom = (struct geometry *) ggg; @@ -692,7 +693,7 @@ check_floors_rad(ldouble *pp, int whichvel,void *ggg) #endif //EVOLVEPHOTONNUMBER #endif //RADIATION - +#endif //SKIPALLFLOORS return ret;