From c447718eebbc7d469e75ca5fda9d2afe2fa158ed Mon Sep 17 00:00:00 2001 From: Andrew Chael Date: Thu, 3 Mar 2022 02:33:08 -0500 Subject: [PATCH] added force free velocities as additional primitives --- Makefile | 4 +- PROBLEMS/MONOPOLE_2D/bc.c | 64 +- PROBLEMS/MONOPOLE_2D/define.h | 126 +- PROBLEMS/MONOPOLE_2D/init.c | 102 +- ana.c | 20 + avg.c | 6 + choices.h | 5 + fileop.c | 182 ++- finite.c | 80 +- frames.c | 57 +- ko.c | 5 + ko.h | 22 +- magn.c | 26 +- metric.c | 3 - misc.c | 62 +- mnemonics.h | 7 + mpi.c | 72 + p2u.c | 125 +- physics.c | 247 ++-- postproc.c | 76 +- problem.c | 6 +- problem.h | 20 +- rad.c | 96 +- relele.c | 74 +- silo.c | 2559 ++++++++++++++++----------------- u2p.c | 661 ++++++--- u2prad.c | 26 +- 27 files changed, 2613 insertions(+), 2120 deletions(-) diff --git a/Makefile b/Makefile index a5a37e3..d5b54a6 100644 --- a/Makefile +++ b/Makefile @@ -2,7 +2,7 @@ ifneq ($(SERIAL),1) CC=mpicc -CFLAGS=-O3 -DMPI +CFLAGS=-O3 -DMPI -fcommon else //CC=gcc @@ -13,7 +13,7 @@ else -fsanitize=address -fno-omit-frame-pointer CC=/usr/bin/h5cc -CFLAGS = -O3 -Wno-unused-result -I/usr/lib/gcc/x86_64-linux-gnu/5.4.0/include -I/usr/include/hdf5/serial -Wunused-function -fopenmp -fcommon +CFLAGS = -O3 -Wno-unused-result -I/usr/lib/gcc/x86_64-linux-gnu/5.4.0/include -I/usr/include/hdf5/serial -Wunused-function -w -fopenmp -fcommon endif diff --git a/PROBLEMS/MONOPOLE_2D/bc.c b/PROBLEMS/MONOPOLE_2D/bc.c index a15eaf2..01e15c6 100755 --- a/PROBLEMS/MONOPOLE_2D/bc.c +++ b/PROBLEMS/MONOPOLE_2D/bc.c @@ -2,6 +2,8 @@ //int calc_bc(int ix,int iy,int iz,ldouble t, // ldouble *uu,ldouble *pp,int ifinit,int BCtype) +#define BCCOORDS KSCOORDS + /**********************/ //geometries ldouble gdet_src,gdet_bc,Fx,Fy,Fz,ppback[NV]; @@ -11,16 +13,16 @@ struct geometry geom; fill_geometry(ix,iy,iz,&geom); struct geometry geomBL; -fill_geometry_arb(ix,iy,iz,&geomBL,BLCOORDS); +fill_geometry_arb(ix,iy,iz,&geomBL,BCCOORDS); -ldouble gg[4][5],GG[4][5],ggsrc[4][5],eup[4][4],elo[4][4]; +ldouble gg[4][5],GG[4][5]; pick_g(ix,iy,iz,gg); pick_G(ix,iy,iz,GG); /**********************/ //radius -if(BCtype==XBCHI) //outflow in magn, atm in rad., atm. in HD + if(BCtype==XBCHI) //outflow in magn, atm in rad., atm. in HD { iix=NX-1; iiy=iy; @@ -32,54 +34,62 @@ if(BCtype==XBCHI) //outflow in magn, atm in rad., atm. in HD pp[iv]=get_u(p,iv,iix,iiy,iiz); } + //ANDREW TODO what to do about outer BC + /* //!! begin rescale //first transform to BL - trans_pmhd_coco(pp, pp, MYCOORDS,BLCOORDS, geom.xxvec,&geom,&geomBL); - + trans_pmhd_coco(pp, pp, MYCOORDS,BCCOORDS, geom.xxvec,&geom,&geomBL); + + // boundary cell and scale factors struct geometry geombdBL; - fill_geometry_arb(iix,iiy,iiz,&geombdBL,BLCOORDS); + fill_geometry_arb(iix,iiy,iiz,&geombdBL,BCCOORDS); ldouble rghost = geomBL.xx; ldouble rbound = geombdBL.xx; ldouble deltar = rbound - rghost; - ldouble scale1 = rbound*rbound/rghost/rghost; - ldouble scale2 = rbound/rghost; - - pp[RHO]*=scale1; - pp[UU] *=scale1; + ldouble scale1 = rbound/rghost; //r^-1 + ldouble scale2 = rbound*rbound/rghost/rghost; //r^-2 + ldouble scale3 = rbound*rbound*rbound/rghost/rghost/rghost; //r^-3 -- inital atm + + // bsq = br*br + r^2bth^2 + r^2*sin^2(th)*bph^2 + + pp[RHO]*=scale3;//scale2; + pp[UU] *=scale3;//scale2; #ifdef MAGNFIELD - pp[B1] *=scale1; - pp[B2] *=scale2; - pp[B3] *=scale2; + pp[B1] *=scale2; + pp[B2] *=scale1; + pp[B3] *=scale1; #endif - //pp[VX] *=1.; + pp[VX] *=1.; pp[VY] *=scale1; pp[VZ] *=scale1; #ifdef RADIATION - pp[EE0] *=scale1; - //pp[FX0] *=1.; - pp[FY0] *=scale1; - pp[FZ0] *=scale1; + pp[EE0] *=scale2; + pp[FX0] *=1.; + pp[FY0] *=scale2; + pp[FZ0] *=scale2; #endif //transform back after rescaling - trans_pmhd_coco(pp, pp,BLCOORDS, MYCOORDS, geomBL.xxvec,&geomBL, &geom); - //!! end rescale + trans_pmhd_coco(pp, pp,BCCOORDS, MYCOORDS, geomBL.xxvec,&geomBL, &geom); - //checking for the gas inflow + //!! end rescale + /* + //check for gas inflow ldouble ucon[4]={0.,pp[VX],pp[VY],pp[VZ]}; conv_vels(ucon,ucon,VELPRIM,VEL4,geom.gg,geom.GG); - trans2_coco(geom.xxvec,ucon,ucon,MYCOORDS,BLCOORDS); - if(ucon[1]<0.) //inflow, resseting to atmosphere + trans2_coco(geom.xxvec,ucon,ucon,MYCOORDS,BCCOORDS); + if(ucon[1]<0.) //inflow, reset to zero velocity { //atmosphere in rho,uint and velocities and zero magn. field //set_hdatmosphere(pp,xxvec,gg,GG,4); ucon[1]=0.; - trans2_coco(geomBL.xxvec,ucon,ucon,BLCOORDS,MYCOORDS); + trans2_coco(geomBL.xxvec,ucon,ucon,BCCOORDS,MYCOORDS); conv_vels(ucon,ucon,VEL4,VELPRIM,geom.gg,geom.GG); pp[VX]=ucon[1]; pp[VY]=ucon[2]; pp[VZ]=ucon[3];//atmosphere in rho,uint and velocities and zero magn. field } - + */ + p2u(pp,uu,&geom); return 0; } @@ -91,7 +101,7 @@ if(BCtype==XBCHI) //outflow in magn, atm in rad., atm. in HD for(iv=0;iv6 ANDREW -#define RMAX 1000. +#define RMAX 100. #define MKS1R0 MKSR0 -#ifdef myMKS2COORDS //modified Kerr-Shild +#ifdef myMKS1COORDS //modified Kerr-Shild +#define MYCOORDS MKS1COORDS +#define MINX (log(RMIN - MKSR0)) +#define MAXX (log(RMAX - MKSR0)) +#define MINY 0.01*M_PI/2. +#define MAXY M_PI - 0.01*M_PI/2. +#define TNX 128//256 +#define TNY 128//256 +#define TNZ 1 + +#elif defined(myMKS2COORDS) //modified Kerr-Shild #define MYCOORDS MKS2COORDS +#define MKSR0 -4.*RHOR +#define MKSH0 0.8 + #define MINX (log(RMIN - MKSR0)) #define MAXX (log(RMAX - MKSR0)) -//#define MINY (0.1*Pi/2.) -//#define MAXY (Pi-0.1*Pi/2.) #define MINY (0.001) #define MAXY (1.-0.001) -#define TNX 512 -#define TNY 40 +#define TNX 128 +#define TNY 128 #define TNZ 1 #else //Schwarzschild @@ -95,26 +125,26 @@ #define PHIWEDGE (M_PI/2.) #define MINZ (-PHIWEDGE/2.) #define MAXZ (PHIWEDGE/2.) -//#define MINZ -1. -//#define MAXZ 1. #define SPECIFIC_BC -#define PEROIDIC_ZBC +#define PERIODIC_ZBC //number of tiles #define NTX 4//48 -#define NTY 1//24 +#define NTY 2//24 #define NTZ 1 /************************************/ //output /************************************/ -#define OUTCOORDS KERRCOORDS +#define OUTCOORDS KSCOORDS //MKS1COORDS #define OUTVEL VEL4 #define ALLSTEPSOUTPUT 0 #define NSTEPSTOP 1.e10 -#define NOUTSTOP 5000 +#define NOUTSTOP 100 #define SILOOUTPUT 1 +#define PRIMOUTPUT 1 +//#define PRIMOUTPUTINMYCOORDS #define SCAOUTPUT 0 #define SILO2D_XZPLANE @@ -122,15 +152,7 @@ //common physics / atmosphere /************************************/ #define GAMMA (4./3.) -#define NODONUT 0 -#define INFLOWING 0 -#define ELL 3.85 -#define URIN (0.) -#define KKK 7127. -#define UTPOT .965 -#define DTOUT1 1.e0 +#define DTOUT1 1 #define DTOUT2 1.e0 #define RHOATMMIN 1.e-20 #define UINTATMMIN 1.e-20 -//#define UINTATMMIN (calc_PEQ_ufromTrho(1.e10,RHOATMMIN)) -//#define ERADATMMIN (calc_LTE_EfromT(3.e6)/10) diff --git a/PROBLEMS/MONOPOLE_2D/init.c b/PROBLEMS/MONOPOLE_2D/init.c index 4436176..7768ee1 100755 --- a/PROBLEMS/MONOPOLE_2D/init.c +++ b/PROBLEMS/MONOPOLE_2D/init.c @@ -1,86 +1,84 @@ ldouble rho,mx,my,mz,m,E,uint,Fx,Fy,Fz,pLTE; -ldouble xx,yy,zz; -ldouble uu[NV],xxvec[4],xxvecBL[4]; -ldouble pp[NV],ppback[NV],T; +ldouble uu[NV],pp[NV],ppback[NV],T; //geometries -get_xx(ix,iy,iz,xxvec); -coco_N(xxvec,xxvecBL,MYCOORDS,BLCOORDS); -xx=xxvecBL[1]; -yy=xxvecBL[2]; -zz=xxvecBL[3]; - -ldouble gg[4][5],GG[4][5],eup[4][4],elo[4][4],tlo[4][4]; -pick_g(ix,iy,iz,gg); -pick_G(ix,iy,iz,GG); - struct geometry geom; fill_geometry(ix,iy,iz,&geom); struct geometry geomBL; -fill_geometry_arb(ix,iy,iz,&geomBL,KERRCOORDS); - -ldouble ggBL[4][5],GGBL[4][5]; -calc_g_arb(xxvecBL,ggBL,KERRCOORDS); -calc_G_arb(xxvecBL,GGBL,KERRCOORDS); - -//ldouble eupBL[4][4],eloBL[4][4]; -//ldouble tupBL[4][4],tloBL[4][4]; -//calc_tetrades(ggBL,tupBL,tloBL,KERRCOORDS); -//calc_ZAMOes(ggBL,eupBL,eloBL,KERRCOORDS); +fill_geometry_arb(ix,iy,iz,&geomBL,BLCOORDS); -ldouble r=geomBL.xx; -ldouble th=geomBL.yy; +ldouble r = geomBL.xx; +ldouble th = geomBL.yy; ldouble rhor = (1. + sqrt(1. - BHSPIN*BHSPIN)) ; -rho = RHOATMMIN + (r/10./rhor)/pow(r,4)/B2RHORATIOMAX; -uint = UINTATMMIN + (r/10./rhor)/pow(r,4)/B2RHORATIOMAX; - -//4-velocity in BL transformed to MYCOORDS -//ldouble ucon[4]={0.,0.,0.,0.}; -//ucon[0] = -1/sqrt(-GGBL[0][0]) -//ucon[3] = GG[0][3]*ucon[0] -//conv_vels(ucon,ucon,VEL3,VELPRIM,ggBL,GGBL); +rho = RHOATMMIN + (r/10./rhor)/pow(r,4)/B2RHORATIOMAXINIT; +uint = UINTATMMIN + (r/10./rhor)/pow(r,4)/B2UURATIOMAXINIT; pp[0]=rho; pp[1]=uint; -pp[2]=0.; //zero in co-rotating frame +pp[2]=0.; //zero initial velocity in co-rotating frame pp[3]=0.; pp[4]=0.; -#ifdef MAGNFIELD//setting them zero not to break the following coordinate transformation - pp[B1]=pp[B2]=pp[B3]=0.; - pp[B1]=1./(r*r); +//entropy +pp[5]=calc_Sfromu(pp[0],pp[1],geom.ix,geom.iy,geom.iz); + +#ifdef MAGNFIELD // set to zero at first +pp[B1]=pp[B2]=pp[B3]=0.; #endif - //transforming primitives from BL to MYCOORDS - trans_pall_coco(pp, pp, KERRCOORDS, MYCOORDS,xxvecBL,&geomBL,&geom); +//transform primitives from BL to MYCOORDS +//trans_pall_coco(pp, pp, BLCOORDS, MYCOORDS,xxvecBL,&geomBL,&geom); -#ifdef MAGNFIELD //MYCOORDS vector potential to calculate B's - ldouble Acov[4]; - Acov[0]=Acov[1]=Acov[2]=0.; - Acov[3]=(1.-cos(th));///(r*sin(th)); - -//pp[B1]=0.; -//pp[B2]=0.; -//pp[B3]=Acov[3]; +#ifdef MAGNFIELD //calculate vector potential +ldouble r_mag,th_mag; + +#ifdef INIT_MAGN_CORNERS +// ANDREW define vector potential at corners not cell centers +// TODO: Issues going to grid edge! We assume we are well inside grid! +ldouble xxvec_c[4],xxvecBL_c[4]; +xxvec_c[0] = global_time; +xxvec_c[1] = get_xb(ix,0); +xxvec_c[2] = get_xb(iy,1); +xxvec_c[3] = get_xb(iz,2); +coco_N(xxvec_c,xxvecBL_c,MYCOORDS,BLCOORDS); + +ldouble r_c=xxvecBL_c[1]; +ldouble th_c=xxvecBL_c[2]; + +r_mag=r_c; +th_mag=th_c; +#else +r_mag=r; +th_mag=th; #endif -//entropy -pp[5]=calc_Sfromu(pp[0],pp[1],geom.ix,geom.iy,geom.iz); +ldouble Acov[4]; +Acov[0]=Acov[1]=Acov[2]=0.; +Acov[3]=(1.-cos(th_mag)); -//to conserved +// Vector potential A is temporarily saved in pp[B1], pp[B2], pp[B3]. +// These will be replaced by B components immediately +pp[B1]=Acov[1]; +pp[B2]=Acov[2]; +pp[B3]=Acov[3]; +#endif + +//all primitives to conserved p2u(pp,uu,&geom); +/***********************************************/ +// set the final values /***********************************************/ int iv; for(iv=0;iv %e %e %e\n",ix,r,th,ph); - print_4vector(xxvec); - print_4vector(xxvecBL); - print_metric(geomBL.gg); - print_metric(geom.gg); - print_4vector(ucon); - */ - - //conv_vels(ucon,ucon,VELPRIM,VEL4,geomBL.gg,geomBL.GG); - //trans2_coco(geomBL.xxvec,ucon,ucon,BLCOORDS,MYCOORDS); - //conv_vels(ucon,ucon,VEL4,VELPRIM,geom.gg,geom.GG); - // print_4vector(ucon); - //getch(); pp[VX]=ucon[1]; pp[VY]=ucon[2]; @@ -6111,12 +6094,7 @@ correct_polaraxis_3d() vr/=sqrt(geom.gg[1][1]); vth/=sqrt(geom.gg[2][2]); vph/=sqrt(geom.gg[3][3]); - //if(ic==0 && ix==10) - // printf("1 %d > %e %e | %e %e %e\n",iy,vth,sqrt(geom.gg[2][2]),vx,vy,vz); ucon[1]=vr; ucon[2]=vth; ucon[3]=vph; - //conv_vels(ucon,ucon,VELPRIM,VEL4,geomBL.gg,geomBL.GG); - //trans2_coco(geomBL.xxvec,ucon,ucon,BLCOORDS,MYCOORDS); - //conv_vels(ucon,ucon,VEL4,VELPRIMRAD,geom.gg,geom.GG); pp[FX]=ucon[1]; pp[FY]=ucon[2]; diff --git a/frames.c b/frames.c index 2121d58..d7fdb3b 100644 --- a/frames.c +++ b/frames.c @@ -180,9 +180,9 @@ calc_Lorentz_lab2ff(ldouble *pp, ldouble gg[][5], ldouble GG[][5], ldouble L[][4 //calculating the four-velocity of fluid in lab frame ldouble utcon[4],ucon[4],ucov[4],vpr[3]; - utcon[1]=pp[2]; - utcon[2]=pp[3]; - utcon[3]=pp[4]; + utcon[1]=pp[VX]; + utcon[2]=pp[VY]; + utcon[3]=pp[VZ]; conv_vels_both(utcon,ucon,ucov,VELPRIM,VEL4,gg,GG); if(verbose>0) print_4vector(ucon); @@ -279,9 +279,9 @@ calc_Lorentz_ff2lab(ldouble *pp, ldouble gg[][5], ldouble GG[][5], ldouble L[][4 //calculating the four-velocity of fluid in lab frame ldouble wcon[4],wtcon[4],wcov[4]; - wtcon[1]=pp[2]; - wtcon[2]=pp[3]; - wtcon[3]=pp[4]; + wtcon[1]=pp[VX]; + wtcon[2]=pp[VY]; + wtcon[3]=pp[VZ]; conv_vels_both(wtcon,wcon,wcov,VELPRIM,VEL4,gg,GG); if(verbose>0) print_4vector(wcon); @@ -468,9 +468,9 @@ boost22_rf2lab(ldouble T1[][4],ldouble T2[][4],ldouble *pp0,ldouble gg[][5],ldou //artificial and temporary substitution ldouble urf[4]={0.,pp[FX],pp[FY],pp[FZ]}; conv_vels(urf,urf,VELPRIMRAD,VELPRIM,gg,GG); - pp[2]=urf[1]; - pp[3]=urf[2]; - pp[4]=urf[3]; + pp[VX]=urf[1]; + pp[VY]=urf[2]; + pp[VZ]=urf[3]; //Lorentz transformation matrix ldouble L[4][4]; @@ -547,9 +547,9 @@ boost22_lab2rf(ldouble T1[][4],ldouble T2[][4],ldouble *pp0,ldouble gg[][5],ldou //artificial and temporary substitution ldouble urf[4]={0.,pp[FX],pp[FY],pp[FZ]}; conv_vels(urf,urf,VELPRIMRAD,VELPRIM,gg,GG); - pp[2]=urf[1]; - pp[3]=urf[2]; - pp[4]=urf[3]; + pp[VX]=urf[1]; + pp[VY]=urf[2]; + pp[VZ]=urf[3]; //Lorentz transformation matrix ldouble L[4][4]; @@ -719,9 +719,9 @@ boost2_lab2rf(ldouble A1[4],ldouble A2[4],ldouble *pp0,ldouble gg[][5],ldouble G //artificial and temporary substitution ldouble urf[4]={0.,pp[FX],pp[FY],pp[FZ]}; conv_vels(urf,urf,VELPRIMRAD,VELPRIM,gg,GG); - pp[2]=urf[1]; - pp[3]=urf[2]; - pp[4]=urf[3]; + pp[VX]=urf[1]; + pp[VY]=urf[2]; + pp[VZ]=urf[3]; if(verbose>0) print_4vector(A1); @@ -889,9 +889,6 @@ multiply2(ldouble *u1,ldouble *u2,ldouble A[][4]) for(i=0;i<4;i++) ut[i]=u1[i]; -#ifdef APPLY_OMP_SIMD - //#pragma omp simd -#endif for(i=0;i<4;i++) { u2[i]=0.; @@ -1456,9 +1453,6 @@ indices_2221(ldouble T1[][4],ldouble T2[][4],ldouble gg[][5]) for(i=0;i<4;i++) { int j; -#ifdef APPLY_OMP_SIMD - //#pragma omp simd -#endif for(j=0;j<4;j++) { Tt[i][j]=0.; @@ -1471,9 +1465,6 @@ indices_2221(ldouble T1[][4],ldouble T2[][4],ldouble gg[][5]) for(j=0;j<4;j++) { int k; -#ifdef APPLY_OMP_SIMD - //#pragma omp simd -#endif for(k=0;k<4;k++) { Tt[i][j]+=T1[i][k]*gg[k][j]; @@ -1484,9 +1475,6 @@ indices_2221(ldouble T1[][4],ldouble T2[][4],ldouble gg[][5]) for(i=0;i<4;i++) { int j; -#ifdef APPLY_OMP_SIMD - //#pragma omp simd -#endif for(j=0;j<4;j++) { T2[i][j]=Tt[i][j]; @@ -1505,17 +1493,11 @@ indices_12(ldouble A1[4],ldouble A2[4],ldouble GG[][5]) { int i; -#ifdef APPLY_OMP_SIMD - //#pragma omp simd -#endif for(i = 0; i < 4; i++) { A2[i] = 0.; } -#ifdef APPLY_OMP_SIMD - //#pragma omp simd -#endif for(i = 0; i < 4; i++) { int k; @@ -1540,9 +1522,6 @@ indices_21(ldouble A1[4],ldouble A2[4],ldouble gg[][5]) int i; ldouble At[4]; -#ifdef APPLY_OMP_SIMD - //#pragma omp simd -#endif for(i=0;i<4;i++) { At[i]=0.; @@ -1551,18 +1530,12 @@ indices_21(ldouble A1[4],ldouble A2[4],ldouble gg[][5]) for(i=0;i<4;i++) { int j; -#ifdef APPLY_OMP_SIMD - //#pragma omp simd -#endif for(j = 0; j < 4; j++) { At[i] += A1[j] * gg[i][j]; } } -#ifdef APPLY_OMP_SIMD - //#pragma omp simd -#endif for (i=0; i<4; i++) { A2[i] = At[i]; diff --git a/ko.c b/ko.c index ed89508..3ec8980 100644 --- a/ko.c +++ b/ko.c @@ -112,6 +112,7 @@ main(int argc, char **argv) #endif //print scalings GU->CGS + if(PROCID==0) printf("NVHD %d NVMHD %d NV %d\n",NVHD,NVMHD,NV); if(PROCID==0) print_scalings(); //exit(1); @@ -468,6 +469,10 @@ main(int argc, char **argv) fprint_simplefile(tstart,nfout1,folder,"sim"); #endif + #if(PRIMOUTPUT!=0) // primitive file + fprint_primitive_file(tstart,nfout1,folder,"prim"); + #endif + #if(RELELSPECTRUMOUTPUT==1) //nonthermal spectrum fprint_relel_spectrum(tstart,NTH_SPEC_IX,NTH_SPEC_IY,NTH_SPEC_IZ,nfout1,folder,"spe",0); #endif diff --git a/ko.h b/ko.h index bfad495..b05c011 100644 --- a/ko.h +++ b/ko.h @@ -345,10 +345,6 @@ ldouble *u,*x,*ucent,*xb,*xout, *xbout_xface, *xbout_yface, *xbout_zface, *flbx,*flby,*flbz,*flLx,*flRx,*flLy,*flRy,*flLz,*flRz, *flbx2,*flby2,*flbz2,*flLx2,*flRx2,*flLy2,*flRy2,*flLz2,*flRz2, *gKr; - //*emuup,*emulo,*emuupbx,*emulobx,*emuupby,*emuloby,*emuupbz,*emulobz, - //*emuup2,*emulo2,*emuupbx2,*emulobx2,*emuupby2,*emuloby2,*emuupbz2,*emulobz2, - //*tmuup,*tmulo,*tmuupbx,*tmulobx,*tmuupby,*tmuloby,*tmuupbz,*tmulobz, - //*tmuup2,*tmulo2,*tmuupbx2,*tmulobx2,*tmuupby2,*tmuloby2,*tmuupbz2,*tmulobz2; int *cellflag; @@ -1008,10 +1004,10 @@ int print_u(ldouble *p); // p2u.c ////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////// -int p2u(ldouble *p, ldouble *u, void *ggg); -int p2u_mhd(ldouble *p, ldouble *u, void *ggg); +int p2u(ldouble *pp, ldouble *uu, void *ggg); +int p2u_mhd(ldouble *pp, ldouble *uu, void *ggg); ldouble calc_utp1(ldouble *vcon, ldouble *ucon, void *ggg); -int p2u_mhd_nonrel(ldouble *p, ldouble *u, void *ggg); +int p2u_mhd_nonrel(ldouble *pp, ldouble *uu, void *ggg); int p2u_rad(ldouble *pp,ldouble *uu,void *ggg); int p2avg(int ix,int iy,int iz,ldouble *avg); int test_maginv(); @@ -1028,6 +1024,7 @@ 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 count_entropy(int *n, int *n2); @@ -1086,9 +1083,6 @@ int set_x(int ic, int idim,ldouble val); int set_xb(int ic, int idim,ldouble val); ldouble calc_xb(int i,int idim); int calc_bc(int ix,int iy,int iz,ldouble t, ldouble *uu,ldouble *pp,int ifinit,int BCtype); -//int set_g(ldouble* uarr,int i,int j,int ix,int iy,int iz,ldouble value); -//int set_T(ldouble* uarr,int i,int j,int ix,int iy,int iz,ldouble value); -//int set_Tfull(ldouble* uarr,int i,int j,int ix,int iy,int iz,ldouble value); int set_ub(ldouble* uarr,int iv,int ix,int iy,int iz,ldouble value,int idim); int set_gb(ldouble* uarr,int i,int j,int ix,int iy,int iz,ldouble value,int idim); int set_Tb(ldouble* uarr,int i,int j,int ix,int iy,int iz,ldouble value,int idim); @@ -1144,10 +1138,10 @@ int mix_entropies(ldouble dt); // magn.c ///////////////////////////////////////////////////// /////////////////////////////////////////////////////////////// -void calc_bcon_bcov_bsq_from_4vel(ldouble *pr, ldouble *ucon, ldouble *ucov, void* ggg, +void calc_bcon_bcov_bsq_from_4vel(ldouble *pp, ldouble *ucon, ldouble *ucov, void* ggg, ldouble *bcon, ldouble *bcov, double *bsq); -void calc_bcon_4vel(double *pr, double *ucon, double *ucov, double *bcon); -void calc_Bcon_4vel(double *pr, double *ucon, double *bcon, double *Bcon); +void calc_bcon_4vel(double *pp, double *ucon, double *ucov, double *bcon); +void calc_Bcon_4vel(double *pp, double *ucon, double *bcon, double *Bcon); void calc_bcon_prim(double *pp, double *bcon, void* ggg); void calc_Bcon_prim(double *pp, double *bcon,double *Bcon, void* ggg); @@ -1306,6 +1300,7 @@ int f_general_source_term(int ix, int iy, int iz,ldouble *ss); int f_flux_prime(ldouble *pp, int idim, int ix, int iy, int iz,ldouble *ff,int lr); int calc_Tij(ldouble *pp, void* ggg, ldouble T[][4]); +int calc_Tij_ff(ldouble *pp, void* ggg, ldouble T[][4]); ldouble calc_ufromS(ldouble S,ldouble rho,int ix,int iy,int iz); ldouble calc_TfromS(ldouble S,ldouble rho,int ix,int iy,int iz); @@ -1461,6 +1456,7 @@ int fprint_simplecart(ldouble t, int nfile, char* folder,char* prefix); int fprint_simplesph(ldouble t, int nfile, char* folder,char* prefix); int fprint_simple_phiavg(ldouble t, int nfile, char* folder,char* prefix); int fprint_simple_phicorr(ldouble t, int nfile, char* folder,char* prefix); +int fprint_primitive_file(ldouble t, int nfile, char* folder, char* prefix); int fprint_restartfile_mpi_hdf5(ldouble t, char* folder); int fprint_restartfile_serial_hdf5(ldouble t, char* folder); diff --git a/magn.c b/magn.c index 4387171..e348db8 100644 --- a/magn.c +++ b/magn.c @@ -8,7 +8,7 @@ /* calculate both magnetic field four-vectors and bsq knowing gas four-velocity ucov */ //*********************************************************************** -void calc_bcon_bcov_bsq_from_4vel(ldouble *pr, ldouble *ucon, ldouble *ucov, void* ggg, +void calc_bcon_bcov_bsq_from_4vel(ldouble *pp, ldouble *ucon, ldouble *ucov, void* ggg, ldouble *bcon, ldouble *bcov, ldouble *bsq) { @@ -17,23 +17,20 @@ void calc_bcon_bcov_bsq_from_4vel(ldouble *pr, ldouble *ucon, ldouble *ucov, voi = (struct geometry *) ggg; // First calculate bcon0 - bcon[0] = pr[B1]*ucov[1] + pr[B2] * ucov[2] + pr[B3] * ucov[3] ; + bcon[0] = pp[B1]*ucov[1] + pp[B2] * ucov[2] + pp[B3] * ucov[3] ; // Then spatial components of bcon #ifdef NONRELMHD for(j = 1; j < 4; j++) - bcon[j] = pr[B1-1+j]; //b^i=B^i + bcon[j] = pp[B1-1+j]; //b^i=B^i #else // relativistic case ldouble u0inv = 1. / ucon[0]; -#ifdef APPLY_OMP_SIMD - //#pragma omp simd -#endif for(j=1;j<4;j++) - bcon[j] = (pr[B1-1+j] + bcon[0] * ucon[j]) * u0inv ; + bcon[j] = (pp[B1-1+j] + bcon[0] * ucon[j]) * u0inv ; #endif //NONRELMHD @@ -48,25 +45,22 @@ void calc_bcon_bcov_bsq_from_4vel(ldouble *pr, ldouble *ucon, ldouble *ucov, voi /* calculate magnetic field four-vector knowing gas four-velocity ucov */ //*********************************************************************** -void calc_bcon_4vel(double *pr, double *ucon, double *ucov, double *bcon) +void calc_bcon_4vel(double *pp, double *ucon, double *ucov, double *bcon) { int j; - bcon[0] = pr[B1]*ucov[1] + pr[B2]*ucov[2] + pr[B3]*ucov[3] ; + bcon[0] = pp[B1]*ucov[1] + pp[B2]*ucov[2] + pp[B3]*ucov[3] ; #ifdef NONRELMHD for(j=1;j<4;j++) - bcon[j] = pr[B1-1+j]; //b^i=B^i + bcon[j] = pp[B1-1+j]; //b^i=B^i #else // relativistic case ldouble u0inv = 1. / ucon[0]; -#ifdef APPLY_OMP_SIMD - //#pragma omp simd -#endif for(j=1;j<4;j++) - bcon[j] = (pr[B1-1+j] + bcon[0]*ucon[j]) * u0inv ; + bcon[j] = (pp[B1-1+j] + bcon[0]*ucon[j]) * u0inv ; #endif //NONRELMHD @@ -77,8 +71,8 @@ void calc_bcon_4vel(double *pr, double *ucon, double *ucov, double *bcon) //*********************************************************************** /* calculate B^i from bcon and gas four-velocity ucon */ //*********************************************************************** - -void calc_Bcon_4vel(double *pr, double *ucon, double *bcon, double *Bcon) +// ANDREW only used once +void calc_Bcon_4vel(double *pp, double *ucon, double *bcon, double *Bcon) { int j; diff --git a/metric.c b/metric.c index 5b20668..a56a2d2 100644 --- a/metric.c +++ b/metric.c @@ -3701,9 +3701,6 @@ dxdx_MKS22KS(ldouble *xx, ldouble dxdx[][4]) #endif int i; -#ifdef APPLY_OMP_SIMD - //#pragma omp simd -#endif for(i=0;i<4;i++) { int j; diff --git a/misc.c b/misc.c index 9df0712..79f5c96 100644 --- a/misc.c +++ b/misc.c @@ -87,15 +87,7 @@ void initialize_constants() two_third = 2. / 3.; log_2p6 = log(2.6); one_over_log_2p6 = 1. / log_2p6; - - /* - printf("Testing gammainterp\n"); - printf("%e %e %e \n",calc_meanlorentz(1.e-3),calc_meanlorentz(1.1e-3),calc_meanlorentz(5.e-3)); - printf("%e %e %e \n",calc_meanlorentz(0.74),calc_meanlorentz(1.23),calc_meanlorentz(7.54)); - printf("%e %e %e \n",calc_meanlorentz(14.3),calc_meanlorentz(77.3),calc_meanlorentz(144.)); - printf("%e %e %e \n",calc_meanlorentz(555.),calc_meanlorentz(999.),calc_meanlorentz(1000.)); - exit(-1); - */ + // Coordinate specific factors #if (MYCOORDS==JETCOORDS) //printf("Finding hypx1out\n"); @@ -107,10 +99,8 @@ void initialize_constants() #ifdef CYLINDRIFY set_cyl_params(); #endif - //ANDREW -- diagnostics for new jet coordinates/metric - /* //printf("%.7f %.7f %.7f\n",hypx1in,hypx1brk,hypx1out); @@ -189,6 +179,16 @@ void initialize_constants() exit(-1); */ + + /* + printf("Testing gammainterp\n"); + printf("%e %e %e \n",calc_meanlorentz(1.e-3),calc_meanlorentz(1.1e-3),calc_meanlorentz(5.e-3)); + printf("%e %e %e \n",calc_meanlorentz(0.74),calc_meanlorentz(1.23),calc_meanlorentz(7.54)); + printf("%e %e %e \n",calc_meanlorentz(14.3),calc_meanlorentz(77.3),calc_meanlorentz(144.)); + printf("%e %e %e \n",calc_meanlorentz(555.),calc_meanlorentz(999.),calc_meanlorentz(1000.)); + exit(-1); + */ + #endif return; @@ -356,6 +356,34 @@ if (NTZ % 2 != 0) printf("\n!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n\n"); } #endif + +#ifdef FORCEFREE + printf("WARNING! FORCEFREE is untested and does not account for current sheets! use only with monopole!\n"); + #ifdef NONRELMHD + printf("FORCEFREE does not work with NONRELMHD!\n"); + exit(-1); + #endif + #ifdef RADIATION + printf("FORCEFREE does not work yet with RADIATION!\n"); + exit(-1); + #endif + #ifndef MAGNFIELD + printf("FORCEFREE requires MAGNFIELD!\n"); + exit(-1); + #endif + #if (VELPRIM!=VELR) + printf("FORCEFREE requries VELPRIM==VELR!\n"); + exit(-1); + #endif + #if(GDETIN==0) + printf("FORCEFREE does not work with GDETIN==0!\n"); + exit(-1); + #endif + #ifdef CORRECT_POLARAXIS_3D + printf("FORCEFREE does not work with CORRECT_POLARAXIS_3D!\n"); + exit(-1); + #endif +#endif #ifdef PWPOTENTIAL printf("PWPOTENTIAL has been removed!\n"); @@ -366,7 +394,7 @@ if (NTZ % 2 != 0) printf("NCOMPTONIZATION has been replaced by EVOLVEPHOTONNUMBER!\n"); exit(-1); #endif - + #ifdef RADIATION #if defined(COMPTONIZATIONFLAG) if (PROCID == 0) @@ -1447,6 +1475,11 @@ print_primitives(ldouble *p) printf("B^1 = %.15e\n",p[B1]); printf("B^2 = %.15e\n",p[B2]); printf("B^3 = %.15e\n",p[B3]); +#ifdef FORCEFREE + printf("u^1 (ff) = %.15e\n",p[VXFF]); + printf("u^2 (ff) = %.15e\n",p[VYFF]); + printf("u^3 (ff) = %.15e\n",p[VZFF]); +#endif #endif #ifdef RADIATION printf("Erf = %.15e\n",p[EE0]); @@ -1490,6 +1523,11 @@ print_conserved(ldouble *u) printf("B^1 = %.15e\n",u[B1]); printf("B^2 = %.15e\n",u[B2]); printf("B^3 = %.15e\n",u[B3]); +#ifdef FORCEFREE + printf("T^t_1 (ff) = %.15e\n",u[VXFF]); + printf("T^t_2 (ff) = %.15e\n",u[VYFF]); + printf("T^t_3 (ff) = %.15e\n",u[VZFF]); +#endif #endif #ifdef RADIATION printf("R^t_t = %.15e\n",u[EE0]); diff --git a/mnemonics.h b/mnemonics.h index 41205c7..e221ff7 100644 --- a/mnemonics.h +++ b/mnemonics.h @@ -34,6 +34,12 @@ #define B2 (NVHD+1) #define B3 (NVHD+2) +#ifdef FORCEFREE +#define VXFF (NVHD+3) +#define VYFF (NVHD+4) +#define VZFF (NVHD+5) +#endif + #define EE (NVMHD) #define FX (NVMHD+1) #define FY (NVMHD+2) @@ -167,6 +173,7 @@ //frames #define ZAMOFRAME 0 #define FFFRAME 1 +#define DRIFTFRAME 2 //avg quantities #define AVGBSQ (NV+0) diff --git a/mpi.c b/mpi.c index 45c2ebf..ff5d52c 100644 --- a/mpi.c +++ b/mpi.c @@ -1993,7 +1993,11 @@ mpi_savedata() { jp = -j + (2*NY + NG - 1); #ifdef MAGNFIELD //need to flip sign of fluxes and velocities across pole, theta and phi + #ifdef FORCEFREE + if(iv==VY || iv==B2 || iv==VYFF || iv==FY0 || iv==VZ || iv==B3 || iv==VZFF || iv==FZ0) + #else if(iv==VY || iv==B2 || iv==FY0 || iv==VZ || iv==B3 || iv==FZ0) + #endif #else if(iv==VY || iv==FY0 || iv==VZ || iv==FZ0) #endif @@ -2034,7 +2038,11 @@ mpi_savedata() { jp = -j - (NG+1); #ifdef MAGNFIELD //need to flip sign of fluxes and velocities across pole, theta and phi + #ifdef FORCEFREE + if(iv==VY || iv==B2 || iv==VYFF || iv==FY0 || iv==VZ || iv==B3 || iv==VZFF || iv==FZ0) + #else if(iv==VY || iv==B2 || iv==FY0 || iv==VZ || iv==B3 || iv==FZ0) + #endif #else if(iv==VY || iv==FY0 || iv==VZ || iv==FZ0) #endif @@ -2097,7 +2105,11 @@ mpi_savedata() { jp = -j + (2*NY + NG - 1); #ifdef MAGNFIELD //need to flip sign of fluxes and velocities across pole, theta and phi + #ifdef FORCEFREE + if(iv==VY || iv==B2 || iv==VYFF || iv==FY0 || iv==VZ || iv==B3 || iv==VZFF || iv==FZ0) + #else if(iv==VY || iv==B2 || iv==FY0 || iv==VZ || iv==B3 || iv==FZ0) + #endif #else if(iv==VY || iv==FY0 || iv==VZ || iv==FZ0) #endif @@ -2138,7 +2150,11 @@ mpi_savedata() { jp = -j - (NG+1); #ifdef MAGNFIELD + #ifdef FORCEFREE + if(iv==VY || iv==B2 || iv==VYFF || iv==FY0 || iv==VZ || iv==B3 || iv==VZFF || iv==FZ0) + #else if(iv==VY || iv==B2 || iv==FY0 || iv==VZ || iv==B3 || iv==FZ0) + #endif #else if(iv==VY || iv==FY0 || iv==VZ || iv==FZ0) #endif @@ -2180,7 +2196,11 @@ mpi_savedata() { jp = -j + (2*NY + NG - 1); #ifdef MAGNFIELD //need to flip sign of fluxes and velocities across pole, theta and phi + #ifdef FORCEFREE + if(iv==VY || iv==B2 || iv==VYFF || iv==FY0 || iv==VZ || iv==B3 || iv==VZFF || iv==FZ0) + #else if(iv==VY || iv==B2 || iv==FY0 || iv==VZ || iv==B3 || iv==FZ0) + #endif #else if(iv==VY || iv==FY0 || iv==VZ || iv==FZ0) #endif @@ -2221,7 +2241,11 @@ mpi_savedata() { jp = -j - (NG+1); #ifdef MAGNFIELD + #ifdef FORCEFREE + if(iv==VY || iv==B2 || iv==VYFF || iv==FY0 || iv==VZ || iv==B3 || iv==VZFF || iv==FZ0) + #else if(iv==VY || iv==B2 || iv==FY0 || iv==VZ || iv==B3 || iv==FZ0) + #endif #else if(iv==VY || iv==FY0 || iv==VZ || iv==FZ0) #endif @@ -2311,7 +2335,11 @@ mpi_savedata() { jp = -j + (2*NY + NG - 1); #ifdef MAGNFIELD //need to flip sign of fluxes and velocities across pole, theta and phi + #ifdef FORCEFREE + if(iv==VY || iv==B2 || iv==VYFF || iv==FY0 || iv==VZ || iv==B3 || iv==VZFF || iv==FZ0) + #else if(iv==VY || iv==B2 || iv==FY0 || iv==VZ || iv==B3 || iv==FZ0) + #endif #else if(iv==VY || iv==FY0 || iv==VZ || iv==FZ0) #endif @@ -2353,7 +2381,11 @@ mpi_savedata() { jp = -j + (2*NY + NG - 1); #ifdef MAGNFIELD //need to flip sign of fluxes and velocities across pole, theta and phi + #ifdef FORCEFREE + if(iv==VY || iv==B2 || iv==VYFF || iv==FY0 || iv==VZ || iv==B3 || iv==VZFF || iv==FZ0) + #else if(iv==VY || iv==B2 || iv==FY0 || iv==VZ || iv==B3 || iv==FZ0) + #endif #else if(iv==VY || iv==FY0 || iv==VZ || iv==FZ0) #endif @@ -2394,7 +2426,11 @@ mpi_savedata() { jp = -j - (NG+1); #ifdef MAGNFIELD + #ifdef FORCEFREE + if(iv==VY || iv==B2 || iv==VYFF || iv==FY0 || iv==VZ || iv==B3 || iv==VZFF || iv==FZ0) + #else if(iv==VY || iv==B2 || iv==FY0 || iv==VZ || iv==B3 || iv==FZ0) + #endif #else if(iv==VY || iv==FY0 || iv==VZ || iv==FZ0) #endif @@ -2438,7 +2474,11 @@ mpi_savedata() { jp = -j - (NG+1); #ifdef MAGNFIELD + #ifdef FORCEFREE + if(iv==VY || iv==B2 || iv==VYFF || iv==FY0 || iv==VZ || iv==B3 || iv==VZFF || iv==FZ0) + #else if(iv==VY || iv==B2 || iv==FY0 || iv==VZ || iv==B3 || iv==FZ0) + #endif #else if(iv==VY || iv==FY0 || iv==VZ || iv==FZ0) #endif @@ -2480,7 +2520,11 @@ mpi_savedata() { jp = -j + (2*NY + NG - 1); #ifdef MAGNFIELD + #ifdef FORCEFREE + if(iv==VY || iv==B2 || iv==VYFF || iv==FY0 || iv==VZ || iv==B3 || iv==VZFF || iv==FZ0) + #else if(iv==VY || iv==B2 || iv==FY0 || iv==VZ || iv==B3 || iv==FZ0) + #endif #else if(iv==VY || iv==FY0 || iv==VZ || iv==FZ0) #endif @@ -2521,7 +2565,11 @@ mpi_savedata() { jp = -j + (2*NY + NG - 1); #ifdef MAGNFIELD + #ifdef FORCEFREE + if(iv==VY || iv==B2 || iv==VYFF || iv==FY0 || iv==VZ || iv==B3 || iv==VZFF || iv==FZ0) + #else if(iv==VY || iv==B2 || iv==FY0 || iv==VZ || iv==B3 || iv==FZ0) + #endif #else if(iv==VY || iv==FY0 || iv==VZ || iv==FZ0) #endif @@ -2563,7 +2611,11 @@ mpi_savedata() { jp = -j - (NG+1); #ifdef MAGNFIELD + #ifdef FORCEFREE + if(iv==VY || iv==B2 || iv==VYFF || iv==FY0 || iv==VZ || iv==B3 || iv==VZFF || iv==FZ0) + #else if(iv==VY || iv==B2 || iv==FY0 || iv==VZ || iv==B3 || iv==FZ0) + #endif #else if(iv==VY || iv==FY0 || iv==VZ || iv==FZ0) #endif @@ -2604,7 +2656,11 @@ mpi_savedata() { jp = -j - (NG+1); #ifdef MAGNFIELD + #ifdef FORCEFREE + if(iv==VY || iv==B2 || iv==VYFF || iv==FY0 || iv==VZ || iv==B3 || iv==VZFF || iv==FZ0) + #else if(iv==VY || iv==B2 || iv==FY0 || iv==VZ || iv==B3 || iv==FZ0) + #endif #else if(iv==VY || iv==FY0 || iv==VZ || iv==FZ0) #endif @@ -2645,7 +2701,11 @@ mpi_savedata() { jp = -j + (2*NY + NG - 1); #ifdef MAGNFIELD + #ifdef FORCEFREE + if(iv==VY || iv==B2 || iv==VYFF || iv==FY0 || iv==VZ || iv==B3 || iv==VZFF || iv==FZ0) + #else if(iv==VY || iv==B2 || iv==FY0 || iv==VZ || iv==B3 || iv==FZ0) + #endif #else if(iv==VY || iv==FY0 || iv==VZ || iv==FZ0) #endif @@ -2686,7 +2746,11 @@ mpi_savedata() { jp = -j + (2*NY + NG - 1); #ifdef MAGNFIELD + #ifdef FORCEFREE + if(iv==VY || iv==B2 || iv==VYFF || iv==FY0 || iv==VZ || iv==B3 || iv==VZFF || iv==FZ0) + #else if(iv==VY || iv==B2 || iv==FY0 || iv==VZ || iv==B3 || iv==FZ0) + #endif #else if(iv==VY || iv==FY0 || iv==VZ || iv==FZ0) #endif @@ -2728,7 +2792,11 @@ mpi_savedata() { jp = -j - (NG+1); #ifdef MAGNFIELD + #ifdef FORCEFREE + if(iv==VY || iv==B2 || iv==VYFF || iv==FY0 || iv==VZ || iv==B3 || iv==VZFF || iv==FZ0) + #else if(iv==VY || iv==B2 || iv==FY0 || iv==VZ || iv==B3 || iv==FZ0) + #endif #else if(iv==VY || iv==FY0 || iv==VZ || iv==FZ0) #endif @@ -2769,7 +2837,11 @@ mpi_savedata() { jp = -j - (NG+1); #ifdef MAGNFIELD + #ifdef FORCEFREE + if(iv==VY || iv==B2 || iv==VYFF || iv==FY0 || iv==VZ || iv==B3 || iv==VZFF || iv==FZ0) + #else if(iv==VY || iv==B2 || iv==FY0 || iv==VZ || iv==B3 || iv==FZ0) + #endif #else if(iv==VY || iv==FY0 || iv==VZ || iv==FZ0) #endif diff --git a/p2u.c b/p2u.c index 9ae7236..1feee58 100644 --- a/p2u.c +++ b/p2u.c @@ -11,12 +11,12 @@ //********************************************************************** int -p2u(ldouble *p, ldouble *u, void *ggg) +p2u(ldouble *pp, ldouble *uu, void *ggg) { - p2u_mhd(p,u,ggg); + p2u_mhd(pp,uu,ggg); #ifdef RADIATION - p2u_rad(p,u,ggg); + p2u_rad(pp,uu,ggg); #endif return 0; @@ -29,11 +29,11 @@ p2u(ldouble *p, ldouble *u, void *ggg) int -p2u_mhd(ldouble *p, ldouble *u, void *ggg) +p2u_mhd(ldouble *pp, ldouble *uu, void *ggg) { #ifdef NONRELMHD - p2u_mhd_nonrel(p,u,ggg); + p2u_mhd_nonrel(pp,uu,ggg); return 0; #endif @@ -50,41 +50,40 @@ p2u_mhd(ldouble *p, ldouble *u, void *ggg) gdetu=1.; #endif - ldouble rho=p[0]; - ldouble uu=p[1]; + ldouble rho=pp[RHO]; + ldouble ugas=pp[UU]; ldouble vcon[4],vcov[4],ucon[4],ucov[4]; ldouble bcon[4]={0.,0.,0.,0.}, bcov[4]={0.,0.,0.,0.}, bsq=0.; vcon[0]=0.; - vcon[1]=p[2]; - vcon[2]=p[3]; - vcon[3]=p[4]; - ldouble S=p[5]; + vcon[1]=pp[VX]; + vcon[2]=pp[VY]; + vcon[3]=pp[VX]; + ldouble S=pp[ENTR]; conv_vels_both(vcon,ucon,ucov,VELPRIM,VEL4,gg,GG); #ifdef MAGNFIELD - calc_bcon_bcov_bsq_from_4vel(p, ucon, ucov, geom, bcon, bcov, &bsq); + calc_bcon_bcov_bsq_from_4vel(pp, ucon, ucov, geom, bcon, bcov, &bsq); #endif //************************************ //hydro part //************************************ - ldouble ut=ucon[0]; + ldouble ut = ucon[0]; ldouble rhout = rho*ut; ldouble Sut = S*ut; - //ANDREW Do we need to ensure gamma is consistent here? ldouble gamma=GAMMA; #ifdef CONSISTENTGAMMA gamma=pick_gammagas(geom->ix,geom->iy,geom->iz); #endif ldouble gammam1=gamma-1.; - ldouble pre=(gamma-1.)*uu; - ldouble w=rho+uu+pre; + ldouble pre=(gamma-1.)*ugas; + ldouble w=rho+ugas+pre; ldouble eta=w+bsq; - ldouble etap = uu+pre+bsq; //eta-rho + ldouble etap=ugas+pre+bsq; //eta-rho ldouble ptot=pre+0.5*bsq; //this computes utp1=1+u_t @@ -96,25 +95,25 @@ p2u_mhd(ldouble *p, ldouble *u, void *ggg) ldouble Ttth =eta*ucon[0]*ucov[2] - bcon[0]*bcov[2]; ldouble Ttph =eta*ucon[0]*ucov[3] - bcon[0]*bcov[3]; - u[0]=gdetu*rhout; - u[1]=gdetu*Tttt; - u[2]=gdetu*Ttr; - u[3]=gdetu*Ttth; - u[4]=gdetu*Ttph; - u[5]=gdetu*Sut; + 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=p[ENTRE]*ut; - u[ENTRE]= gdetu*Seut; - ldouble Siut=p[ENTRI]*ut; - u[ENTRI]= gdetu*Siut; + ldouble Seut=pp[ENTRE]*ut; + uu[ENTRE]= gdetu*Seut; + ldouble Siut=pp[ENTRI]*ut; + uu[ENTRI]= gdetu*Siut; #endif #ifdef RELELECTRONS int ib; for(ib=0;ibix,geom->iy,geom->iz); - //ANDREW -- can we do this consistently instead? - //gamma = 1. + (state->pi+state->pe+state->penth)/(state->ue+state->ui+state->uenth); #endif #endif pgas = (gamma - 1.) * uint; state->gamma = gamma; state->pgas = pgas; - //ANDREW this is in different units than in pp[ENTR]! + //ANDREW gas entropy is in different units than in pp[ENTR]! Sgas = kB_over_mugas_mp*calc_Sfromu(rho, uint, geom->ix,geom->iy,geom->iz); Tgas = calc_PEQ_Teifrompp(pp,&Te,&Ti,geom->ix,geom->iy,geom->iz); @@ -52,15 +50,11 @@ fill_struct_of_state(ldouble *pp, void* ggg, void* sss) ldouble pnth=0.; ldouble unth=0.; #ifdef RELELECTRONS - nnth = calc_relel_ne(pp); - //if(nnth/netot > MAX_RELEL_FRAC_N) nnth = MAX_RELEL_FRAC_N*netot; - + nnth = calc_relel_ne(pp); unth = calc_relel_uint(pp); - //if (unth/uint > MAX_RELEL_FRAC_U) unth = MAX_RELEL_FRAC_U*uint; - pnth = calc_relel_p(pp); - //if (pnth/pgas > MAX_RELEL_FRAC_P) pnth = MAX_RELEL_FRAC_P*pgas; #endif + state->nenth=nnth; state->uenth=unth; state->penth=pnth; @@ -110,7 +104,7 @@ fill_struct_of_state(ldouble *pp, void* ggg, void* sss) #ifdef MAGNFIELD calc_bcon_bcov_bsq_from_4vel(pp, ucon, ucov, geom, bcon, bcov, &bsq); - calc_Bcon_4vel(pp, ucon, bcon, Bcon); + calc_Bcon_4vel(pp,ucon, bcon, Bcon); DLOOPA(i) { @@ -502,7 +496,7 @@ calc_wavespeeds_lr_pure(ldouble *pp,void *ggg,ldouble *aaa) ldouble vtot2; //total characteristic velocity vtot2=cs2 + va2 - cs2*va2; - + #ifdef NONRELMHD //non-rel version ldouble vx,vy,vz,cs,csx,csy,csz; @@ -819,10 +813,11 @@ int f_metric_source_term_arb(ldouble *pp,void *ggg,ldouble *ss) dlgdet[1]=gg[1][4]; //D[gdet,x2]/gdet dlgdet[2]=gg[2][4]; //D[gdet,x3]/gdet - ldouble ut; - ldouble T[4][4]; + + //calculating stress energy tensor components + ldouble T[4][4]; calc_Tij(pp,geom,T); indices_2221(T,T,gg); @@ -840,10 +835,10 @@ int f_metric_source_term_arb(ldouble *pp,void *ggg,ldouble *ss) ldouble rho=pp[RHO]; ldouble u=pp[UU]; ldouble vcon[4],ucon[4]; - vcon[1]=pp[2]; - vcon[2]=pp[3]; - vcon[3]=pp[4]; - ldouble S=pp[5]; + vcon[1]=pp[VX]; + vcon[2]=pp[VY]; + vcon[3]=pp[VZ]; + ldouble S=pp[ENTR]; //converting to 4-velocity conv_vels(vcon,ucon,VELPRIM,VEL4,gg,GG); @@ -862,10 +857,10 @@ int f_metric_source_term_arb(ldouble *pp,void *ggg,ldouble *ss) for(k=0;k<4;k++) for(l=0;l<4;l++) { - ss[1]+=gdetu*T[k][l]*get_gKr(l,0,k,ix,iy,iz); - ss[2]+=gdetu*T[k][l]*get_gKr(l,1,k,ix,iy,iz); - ss[3]+=gdetu*T[k][l]*get_gKr(l,2,k,ix,iy,iz); - ss[4]+=gdetu*T[k][l]*get_gKr(l,3,k,ix,iy,iz); + ss[UU]+=gdetu*T[k][l]*get_gKr(l,0,k,ix,iy,iz); + ss[VX]+=gdetu*T[k][l]*get_gKr(l,1,k,ix,iy,iz); + ss[VY]+=gdetu*T[k][l]*get_gKr(l,2,k,ix,iy,iz); + ss[VZ]+=gdetu*T[k][l]*get_gKr(l,3,k,ix,iy,iz); ss[EE0]+=gdetu*Rij[k][l]*get_gKr(l,0,k,ix,iy,iz); ss[FX0]+=gdetu*Rij[k][l]*get_gKr(l,1,k,ix,iy,iz); ss[FY0]+=gdetu*Rij[k][l]*get_gKr(l,2,k,ix,iy,iz); @@ -884,12 +879,12 @@ int f_metric_source_term_arb(ldouble *pp,void *ggg,ldouble *ss) for(l=1;l<4;l++) { - ss[0]+=-dlgdet[l-1]*rho*ucon[l]; - ss[1]+=-dlgdet[l-1]*(T[l][0]+rho*ucon[l]); - ss[2]+=-dlgdet[l-1]*(T[l][1]); - ss[3]+=-dlgdet[l-1]*(T[l][2]); - ss[4]+=-dlgdet[l-1]*(T[l][3]); - ss[5]+=-dlgdet[l-1]*S*ucon[l]; + ss[RHO]+=-dlgdet[l-1]*rho*ucon[l]; + ss[UU]+=-dlgdet[l-1]*(T[l][0]+rho*ucon[l]); + ss[VX]+=-dlgdet[l-1]*(T[l][1]); + ss[VY]+=-dlgdet[l-1]*(T[l][2]); + ss[VZ]+=-dlgdet[l-1]*(T[l][3]); + ss[ENTR]+=-dlgdet[l-1]*S*ucon[l]; ss[EE0]+=-dlgdet[l-1]*(Rij[l][0]); ss[FX0]+=-dlgdet[l-1]*(Rij[l][1]); ss[FY0]+=-dlgdet[l-1]*(Rij[l][2]); @@ -913,32 +908,59 @@ int f_metric_source_term_arb(ldouble *pp,void *ggg,ldouble *ss) } #endif //GDETIN -#else //ndef RADIATION , pure hydro +#else //no RADIATION , pure hydro //terms with Christoffels for(k=0;k<4;k++) for(l=0;l<4;l++) { - ss[1]+=gdetu*T[k][l]*get_gKr(l,0,k,ix,iy,iz); - ss[2]+=gdetu*T[k][l]*get_gKr(l,1,k,ix,iy,iz); - ss[3]+=gdetu*T[k][l]*get_gKr(l,2,k,ix,iy,iz); - ss[4]+=gdetu*T[k][l]*get_gKr(l,3,k,ix,iy,iz); + ss[UU]+=gdetu*T[k][l]*get_gKr(l,0,k,ix,iy,iz); + ss[VX]+=gdetu*T[k][l]*get_gKr(l,1,k,ix,iy,iz); + ss[VY]+=gdetu*T[k][l]*get_gKr(l,2,k,ix,iy,iz); + ss[VZ]+=gdetu*T[k][l]*get_gKr(l,3,k,ix,iy,iz); } //terms with dloggdet #if (GDETIN==0) for(l=1;l<4;l++) { - ss[0]+=-dlgdet[l-1]*rho*ucon[l]; - ss[1]+=-dlgdet[l-1]*(T[l][0]+rho*ucon[l]); - ss[2]+=-dlgdet[l-1]*(T[l][1]); - ss[3]+=-dlgdet[l-1]*(T[l][2]); - ss[4]+=-dlgdet[l-1]*(T[l][3]); - ss[5]+=-dlgdet[l-1]*S*ucon[l]; + ss[RHO]+=-dlgdet[l-1]*rho*ucon[l]; + ss[UU]+=-dlgdet[l-1]*(T[l][0]+rho*ucon[l]); + ss[VX]+=-dlgdet[l-1]*(T[l][1]); + ss[VY]+=-dlgdet[l-1]*(T[l][2]); + ss[VZ]+=-dlgdet[l-1]*(T[l][3]); + ss[ENTR]+=-dlgdet[l-1]*S*ucon[l]; } #endif #endif //RADIATION +#ifdef FORCEFREE + //calculating stress energy tensor components + ldouble T_ff[4][4]; + calc_Tij_ff(pp,geom,T_ff); + indices_2221(T_ff,T_ff,gg); + + for(ii=0;ii<4;ii++) + for(jj=0;jj<4;jj++) + { + if(isnan(T_ff[ii][jj])) + { + printf("%d %d %e\n",ii,jj,T[ii][jj]); + my_err("nan in force free metric_source_terms\n"); + } + } + + //terms with Christoffels + for(k=0;k<4;k++) + for(l=0;l<4;l++) + { + ss[VXFF]+=gdetu*T_ff[k][l]*get_gKr(l,1,k,ix,iy,iz); + ss[VYFF]+=gdetu*T_ff[k][l]*get_gKr(l,2,k,ix,iy,iz); + ss[VZFF]+=gdetu*T_ff[k][l]*get_gKr(l,3,k,ix,iy,iz); + } + +#endif //FORCEFREE + #ifdef SHEARINGBOX //signs the same despite rho u^i u_t evolved because in koral source terms on lhs @@ -1132,6 +1154,7 @@ int f_general_source_term(int ix, int iy, int iz,ldouble *ss) //*************************************************************** // calculates fluxes at faces //*************************************************************** + int f_flux_prime(ldouble *pp, int idim, int ix, int iy, int iz,ldouble *ff,int lr) { @@ -1153,60 +1176,58 @@ int f_flux_prime(ldouble *pp, int idim, int ix, int iy, int iz,ldouble *ff,int l gdetu=1.; #endif - //calculating Tij + //calculating stress energy tensor Tij ldouble T[4][4]; calc_Tij(pp,&geom,T); indices_2221(T,T,gg);//T^ij --> T^i_j //primitives -#ifdef EVOLVEELECTRONS + #ifdef EVOLVEELECTRONS ldouble Se=pp[ENTRE]; //entropy of electrons ldouble Si=pp[ENTRI]; //entropy of ions -#endif + #endif ldouble rho=pp[RHO]; ldouble u=pp[UU]; - + ldouble S=pp[ENTR]; + ldouble vcon[4],ucon[4],ucov[4],bcon[4],bcov[4],bsq=0.; - - vcon[1]=pp[2]; - vcon[2]=pp[3]; - vcon[3]=pp[4]; - ldouble S=pp[5]; + vcon[1]=pp[VX]; + vcon[2]=pp[VY]; + vcon[3]=pp[VZ]; //converting to 4-velocity conv_vels_both(vcon,ucon,ucov,VELPRIM,VEL4,gg,GG); -#ifdef NONRELMHD + #ifdef NONRELMHD ucon[0]=1.; ucov[0]=-1.; -#endif + #endif -#ifdef MAGNFIELD + #ifdef MAGNFIELD calc_bcon_bcov_bsq_from_4vel(pp, ucon, ucov, &geom, bcon, bcov, &bsq); -#endif + #endif ldouble gamma=GAMMA; #ifdef CONSISTENTGAMMA gamma=pick_gammagas(ix,iy,iz); #endif - ldouble gammam1=gamma-1.; ldouble pre=(gamma-1.)*u; ldouble w=rho+u+pre; ldouble eta=w+bsq; ldouble etap = u+pre+bsq; //eta-rho - int ii, jj, irf; + int ii, jj; for(ii=0;ii<4;ii++) for(jj=0;jj<4;jj++) { if(isnan(T[ii][jj])) { printf("%d > nan tmunu: %d %d %e at %d %d %d\n",PROCID,ii,jj,T[ii][jj],ix+TOI,iy+TOJ,iz+TOK); - printf("%d > nan tmunu: %e %e %e %e\n",PROCID,gamma,pre,w,eta); - print_4vector(ucon); - print_metric(geom.gg); - print_Nvector(pp,NV); + //printf("%d > nan tmunu: %e %e %e %e\n",PROCID,gamma,pre,w,eta); + //print_4vector(ucon); + //print_metric(geom.gg); + //print_Nvector(pp,NV); my_err("nan in flux_prime\n"); exit(1); } @@ -1218,24 +1239,24 @@ int f_flux_prime(ldouble *pp, int idim, int ix, int iy, int iz,ldouble *ff,int l //fluxes //*************************************** //hydro fluxes - ff[0]= gdetu*rho*ucon[idim+1]; - - //ff[1]= gdetu*(T[idim+1][0]+rho*ucon[idim+1]); - //to avoid slow cancellation: - ff[1]= gdetu*(etap*ucon[idim+1]*ucov[0] + rho*ucon[idim+1]*utp1); -#ifdef MAGNFIELD - ff[1]+= -gdetu*bcon[idim+1]*bcov[0]; -#endif - - ff[2]= gdetu*(T[idim+1][1]); - ff[3]= gdetu*(T[idim+1][2]); - ff[4]= gdetu*(T[idim+1][3]); - ff[5]= gdetu*S*ucon[idim+1]; - -#ifdef NONRELMHD - ff[1]= gdetu*T[idim+1][0]; + ff[RHO] = gdetu*rho*ucon[idim+1]; + ff[ENTR]= gdetu*S*ucon[idim+1]; + +#if defined(NONRELMHD) + ff[UU]= gdetu*T[idim+1][0]; +#else + //ff[UU]= gdetu*(T[idim+1][0]+rho*ucon[idim+1]); + //to avoid slow cancellation use the more explicit form: + ff[UU]= gdetu*(etap*ucon[idim+1]*ucov[0] + rho*ucon[idim+1]*utp1); + #ifdef MAGNFIELD + ff[UU]+= -gdetu*bcon[idim+1]*bcov[0]; + #endif #endif + ff[VX]= gdetu*(T[idim+1][1]); + ff[VY]= gdetu*(T[idim+1][2]); + ff[VZ]= gdetu*(T[idim+1][3]); + #ifdef EVOLVEELECTRONS ff[ENTRE]= gdetu*Se*ucon[idim+1]; ff[ENTRI]= gdetu*Si*ucon[idim+1]; @@ -1262,7 +1283,29 @@ int f_flux_prime(ldouble *pp, int idim, int ix, int iy, int iz,ldouble *ff,int l ff[B1]+=suppfac*gdetu*eterm[1]; ff[B2]+=suppfac*gdetu*eterm[2]; ff[B3]+=suppfac*gdetu*eterm[3]; +#endif #endif + +#ifdef FORCEFREE + ldouble T_ff[4][4]; + calc_Tij_ff(pp,&geom,T_ff); + indices_2221(T_ff,T_ff,gg); + + for(ii=0;ii<4;ii++) + for(jj=0;jj<4;jj++) + { + if(isnan(T_ff[ii][jj])) + { + printf("%d > nan tmunu_ff: %d %d %e at %d %d %d\n",PROCID,ii,jj,T_ff[ii][jj],ix+TOI,iy+TOJ,iz+TOK); + my_err("nan in force free flux_prime\n"); + exit(1); + } + } + + ff[VXFF]= gdetu*(T_ff[idim+1][1]); + ff[VYFF]= gdetu*(T_ff[idim+1][2]); + ff[VZFF]= gdetu*(T_ff[idim+1][3]); + #endif #endif @@ -1306,7 +1349,6 @@ calc_Tij(ldouble *pp, void* ggg, ldouble T[][4]) ucov[0]=-1.; #endif - #ifdef MAGNFIELD calc_bcon_bcov_bsq_from_4vel(pp, ucon, ucov, geom, bcon, bcov, &bsq); #else @@ -1318,22 +1360,13 @@ calc_Tij(ldouble *pp, void* ggg, ldouble T[][4]) #ifdef CONSISTENTGAMMA gamma=pick_gammagas(geom->ix,geom->iy,geom->iz); #endif - ldouble gammam1=gamma-1.; ldouble p=(gamma-1.)*uu; ldouble w=rho+uu+p; ldouble eta=w+bsq; ldouble ptot=p+0.5*bsq; -#ifndef NONRELMHD - for(i=0;i<4;i++) -#ifdef APPLY_OMP_SIMD - //#pragma omp simd -#endif - for(j=0;j<4;j++) - T[i][j]=eta*ucon[i]*ucon[j] + ptot*GG[i][j] - bcon[i]*bcon[j]; - -#else +#if defined(NONRELMHD) // non-relativistic ldouble v2=dot3nr(ucon,ucov); @@ -1347,8 +1380,50 @@ calc_Tij(ldouble *pp, void* ggg, ldouble T[][4]) for(i=1;i<4;i++) T[0][i]=T[i][0]=(T[0][0] + ptot) *ucon[i]*ucon[0] + ptot*GG[i][0] - bcon[i]*bcon[0]; -#endif // ifndef NONRELMHD +#else //normal GRMHD + for(i=0;i<4;i++) + for(j=0;j<4;j++) + T[i][j]=eta*ucon[i]*ucon[j] + ptot*GG[i][j] - bcon[i]*bcon[j]; + +#endif + + return 0; +} + +// force-free, or EM part of the stress-energy tensor only +int +calc_Tij_ff(ldouble *pp, void* ggg, ldouble T[][4]) +{ + struct geometry *geom + = (struct geometry *) ggg; + +#ifdef FORCEFREE + ldouble (*gg)[5],(*GG)[5]; + gg=geom->gg; + GG=geom->GG; + + int iv,i,j; + ldouble utcon[4],ucon[4],ucov[4]; + ldouble bcon[4],bcov[4],bsq=0.; + + // ANDREW TODO : which 4-velocity to use here generally? + //get 4-velocity + + utcon[0]=0.; + utcon[1]=pp[VXFF]; // or pp[VX] ? + utcon[2]=pp[VYFF]; // or pp[VY] ? + utcon[3]=pp[VZFF]; // or pp[VZ] ? + + conv_vels_both(utcon,ucon,ucov,VELPRIM,VEL4,gg,GG); + + calc_bcon_bcov_bsq_from_4vel(pp, ucon, ucov, geom, bcon, bcov, &bsq); + + for(i=0;i<4;i++) + for(j=0;j<4;j++) + T[i][j]=bsq*ucon[i]*ucon[j] + 0.5*bsq*GG[i][j] - bcon[i]*bcon[j]; + +#endif // FORCEFREE return 0; } @@ -1694,7 +1769,7 @@ calc_ufromS3rho(ldouble S3,ldouble rho,int type,int ix,int iy,int iz) //********************************************************************** -//updates entropy in the specified cell (p[5]) basing on new primitives +//updates entropy in the specified cell (p[ENTR]) basing on new primitives //or stays with the old one if entropy u2p solver was involved // Ramesh: This function is called by init.c in HIGHZBHS; // we should change all the others diff --git a/postproc.c b/postproc.c index 1bc4d50..4c72aa2 100644 --- a/postproc.c +++ b/postproc.c @@ -416,11 +416,11 @@ int calc_radialprofiles(ldouble profiles[][NX]) } // end if(doingavg) else //on the go from the primitives { - rho=pp[0]; - uint=pp[1]; - utcon[1]=pp[2]; - utcon[2]=pp[3]; - utcon[3]=pp[4]; + rho=pp[RHO]; + uint=pp[UU]; + utcon[1]=pp[VX]; + utcon[2]=pp[VY]; + utcon[3]=pp[VZ]; temp=calc_PEQ_Tfromurho(uint,rho,ix,iy,iz); #ifdef MAGNFIELD @@ -1080,11 +1080,11 @@ int calc_thetaprofiles(ldouble profiles[][NY]) } else { - rho=pp[0]; - uint=pp[1]; - utcon[1]=pp[2]; - utcon[2]=pp[3]; - utcon[3]=pp[4]; + rho=pp[RHO]; + uint=pp[UU]; + utcon[1]=pp[VX]; + utcon[2]=pp[VY]; + utcon[3]=pp[VZ]; #ifdef MAGNFIELD calc_bcon_prim(pp,bcon,&geomBL); @@ -1926,9 +1926,9 @@ calc_local_lum(int ix,int iy,int iz,ldouble *radlum, ldouble *totallum) else { - ucongas[1]=pp[2]; - ucongas[2]=pp[3]; - ucongas[3]=pp[4]; + ucongas[1]=pp[VX]; + ucongas[2]=pp[VY]; + ucongas[3]=pp[VZ]; conv_vels(ucongas,ucongas,VELPRIM,VEL4,geom.gg,geom.GG); rhour = pp[RHO]*ucongas[1]; @@ -2069,9 +2069,9 @@ calc_lum(ldouble radius,int type,ldouble *radlum, ldouble *totallum) else //doingavg { - ucongas[1]=pp[2]; - ucongas[2]=pp[3]; - ucongas[3]=pp[4]; + ucongas[1]=pp[VX]; + ucongas[2]=pp[VY]; + ucongas[3]=pp[VZ]; conv_vels(ucongas,ucongas,VELPRIM,VEL4,geom.gg,geom.GG); rhour = pp[RHO]*ucongas[1]; @@ -2271,9 +2271,9 @@ calc_lum(ldouble radius,int type,ldouble *radlum, ldouble *totallum) //hydro part may be inconsistent! calc_tautot(pp,&geomBL,dxph,tautot); - ucongas[1]=pp[2]; - ucongas[2]=pp[3]; - ucongas[3]=pp[4]; + ucongas[1]=pp[VX]; + ucongas[2]=pp[VY]; + ucongas[3]=pp[VZ]; conv_vels(ucongas,ucongas,VELPRIM,VEL4,geomBL.gg,geomBL.GG); indices_21(ucongas,ucovgas,geomBL.gg); @@ -2456,9 +2456,9 @@ calc_lum_tausurface(ldouble taumax,ldouble *radlum) trans_pall_coco(pp, pp, MYCOORDS,OUTCOORDS, geom.xxvec,&geom,&geomBL); #endif - ucongas[1]=pp[2]; - ucongas[2]=pp[3]; - ucongas[3]=pp[4]; + ucongas[1]=pp[VX]; + ucongas[2]=pp[VY]; + ucongas[3]=pp[VZ]; conv_vels(ucongas,ucongas,VELPRIM,VEL4,geomBL.gg,geomBL.GG); indices_21(ucongas,ucovgas,geomBL.gg); @@ -2917,10 +2917,10 @@ calc_mdot(ldouble radius, int type) pick_g(ix, iy, iz, gg); pick_G(ix, iy, iz, GG); - rho = pp[0]; - ucon[1] = pp[2]; - ucon[2] = pp[3]; - ucon[3] = pp[4]; + rho = pp[RHO]; + ucon[1] = pp[VX]; + ucon[2] = pp[VY]; + ucon[3] = pp[VZ]; conv_vels(ucon, ucon, VELPRIM, VEL4, geom.gg, geom.GG); rhouconr = rho * ucon[1]; @@ -3032,13 +3032,13 @@ calc_lum_proxy(ldouble radius, ldouble theta_min, ldouble theta_max) pick_g(ix, iy, iz, gg); pick_G(ix, iy, iz, GG); - rho = pp[0]; + rho = pp[RHO]; uu = pp[UU]; pressure = uu * (GAMMA - 1.); - ucon[1] = pp[2]; - ucon[2] = pp[3]; - ucon[3] = pp[4]; + ucon[1] = pp[VX]; + ucon[2] = pp[VY]; + ucon[3] = pp[VZ]; conv_vels(ucon, ucon, VELPRIM, VEL4, geom.gg, geom.GG); indices_21(ucon, ucov, geom.gg); @@ -3156,12 +3156,12 @@ calc_Edot(ldouble radius) pick_g(ix, iy, iz, gg); pick_G(ix, iy, iz, GG); - rho = pp[0]; + rho = pp[RHO]; uu = pp[UU]; - ucon[1] = pp[2]; - ucon[2] = pp[3]; - ucon[3] = pp[4]; + ucon[1] = pp[VX]; + ucon[2] = pp[VY]; + ucon[3] = pp[VZ]; conv_vels(ucon, ucon, VELPRIM, VEL4, geom.gg, geom.GG); indices_21(ucon, ucov, geom.gg); @@ -3282,12 +3282,12 @@ calc_Ldot(ldouble radius) pick_g(ix, iy, iz, gg); pick_G(ix, iy, iz, GG); - rho = pp[0]; + rho = pp[RHO]; uu = pp[UU]; - ucon[1] = pp[2]; - ucon[2] = pp[3]; - ucon[3] = pp[4]; + ucon[1] = pp[VX]; + ucon[2] = pp[VY]; + ucon[3] = pp[VZ]; conv_vels(ucon, ucon, VELPRIM, VEL4, geom.gg, geom.GG); indices_21(ucon, ucov, geom.gg); diff --git a/problem.c b/problem.c index b6cbded..4b4f2df 100644 --- a/problem.c +++ b/problem.c @@ -923,7 +923,11 @@ solve_the_problem(ldouble tstart, char* folder) #if(SIMOUTPUT!=0) //sim files fprint_simplefile(tstart,nfout1,folder,"sim"); #endif - + + #if(PRIMOUTPUT!=0) //prim output + fprint_primitive_file(t,nfout1,folder,"prim"); + #endif + #if(RELELSPECTRUMOUTPUT==1) //nonthermal spectrum fprint_relel_spectrum(t,NTH_SPEC_IX,NTH_SPEC_IY,NTH_SPEC_IZ,nfout1,folder,"spe",0); #endif diff --git a/problem.h b/problem.h index bce2c54..ebc2714 100644 --- a/problem.h +++ b/problem.h @@ -144,13 +144,29 @@ //140 RADSURVEY -- radiative parameter survey //141 KEPINF -- INFDISK modified for injecting keplerian //142 PARTIALTDE -- INFDISK modified for partial TDE binding energy distribution non-constant +//143 FFTESTS -- Force Free tests //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 140 +#define PROBLEM 132 //140 + + +#if(PROBLEM==143) +#define PR_DEFINE "PROBLEMS/FFTESTS/define.h" +#define PR_BC "PROBLEMS/FFTESTS/bc.c" +#define PR_INIT "PROBLEMS/FFTESTS/init.c" +#define PR_KAPPAES "PROBLEMS/FFTESTS/kappaes.c" +#define PR_OUT2GIF_2D "PROBLEMS/FFTESTS/out2gif_2d.c" +#define PR_OUT2GIF_1D "PROBLEMS/FFTESTS/out2gif_1d.c" +#define PR_DUMP "PROBLEMS/FFTESTS/dump.c" +#define PR_FINGER "PROBLEMS/FFTESTS/finger.c" +#define PR_TOOLS "PROBLEMS/FFTESTS/tools.c" +#define PR_PREPINIT "PROBLEMS/FFTESTS/prepinit.c" +#define PR_POSTINIT "PROBLEMS/FFTESTS/postinit.c" +#endif #if(PROBLEM==142) @@ -294,7 +310,7 @@ //#define PR_OUT2GIF_1D "PROBLEMS/MONOPOLE_2D/out2gif_1d.c" //#define PR_DUMP "PROBLEMS/MONOPOLE_2D/dump.c" //#define PR_TOOLS "PROBLEMS/MONOPOLE_2D/tools.c" -//#define PR_POSTINIT "PROBLEMS/MAGDONUT/postinit.c" +#define PR_POSTINIT "PROBLEMS/MONOPOLE_2D/postinit.c" #endif #if(PROBLEM==131) diff --git a/rad.c b/rad.c index b5b4538..793bd86 100644 --- a/rad.c +++ b/rad.c @@ -1332,10 +1332,10 @@ implicit_apply_constraints(ldouble *pp, ldouble *uu, ldouble *uu0, void* ggg, in p2u_mhd(pp,uu,geom); //make opposite changes in the RAD conserved quantities - uu[EE0] = uu0[EE0] - (uu[1]-uu0[1]); - uu[FX0] = uu0[FX0] - (uu[2]-uu0[2]); - uu[FY0] = uu0[FY0] - (uu[3]-uu0[3]); - uu[FZ0] = uu0[FZ0] - (uu[4]-uu0[4]); + uu[EE0] = uu0[EE0] - (uu[UU]-uu0[UU]); + uu[FX0] = uu0[FX0] - (uu[VX]-uu0[VX]); + uu[FY0] = uu0[FY0] - (uu[VY]-uu0[VY]); + uu[FZ0] = uu0[FZ0] - (uu[VZ]-uu0[VZ]); //and invert back to primitives u2pret=u2p_rad(uu,pp,geom,corr); @@ -1380,10 +1380,10 @@ implicit_apply_constraints(ldouble *pp, ldouble *uu, ldouble *uu0, void* ggg, in //make opposite changes in the MHD conserved quantities uu[RHO]=uu0[RHO]; - uu[1] = uu0[1] - (uu[EE0]-uu0[EE0]); - uu[2] = uu0[2] - (uu[FX0]-uu0[FX0]); - uu[3] = uu0[3] - (uu[FY0]-uu0[FY0]); - uu[4] = uu0[4] - (uu[FZ0]-uu0[FZ0]); + uu[UU] = uu0[UU] - (uu[EE0]-uu0[EE0]); + uu[VX] = uu0[VX] - (uu[FX0]-uu0[FX0]); + uu[VY] = uu0[VY] - (uu[FY0]-uu0[FY0]); + uu[VZ] = uu0[VZ] - (uu[FZ0]-uu0[FZ0]); //copy B field to uu #ifdef MAGNFIELD @@ -1676,13 +1676,13 @@ f_implicit_lab_4dprim_with_state(ldouble *uu, ldouble *pp, void *sss, //errors in momenta -- always in lab frame if(whichprim==MHD) //mhd-primitives { - f[1] = uu[2] - uu0[2] - dt * gdetu * Gi[1]; - f[2] = uu[3] - uu0[3] - dt * gdetu * Gi[2]; - f[3] = uu[4] - uu0[4] - dt * gdetu * Gi[3]; + f[1] = uu[VX] - uu0[VX] - dt * gdetu * Gi[1]; + f[2] = uu[VY] - uu0[VY] - dt * gdetu * Gi[2]; + f[3] = uu[VZ] - uu0[VZ] - dt * gdetu * Gi[3]; - if(fabs(f[1])>SMALL) err[1]=fabs(f[1])/(1.e-20*uu[UU]+fabs(uu[2])+fabs(uu0[2])+fabs(dt*gdetu*Gi[1])); else err[1]=0.; - if(fabs(f[2])>SMALL) err[2]=fabs(f[2])/(1.e-20*uu[UU]+fabs(uu[3])+fabs(uu0[3])+fabs(dt*gdetu*Gi[2])); else err[2]=0.; - if(fabs(f[3])>SMALL) err[3]=fabs(f[3])/(1.e-20*uu[UU]+fabs(uu[4])+fabs(uu0[4])+fabs(dt*gdetu*Gi[3])); else err[3]=0.; + if(fabs(f[1])>SMALL) err[1]=fabs(f[1])/(1.e-20*uu[UU]+fabs(uu[VX])+fabs(uu0[VX])+fabs(dt*gdetu*Gi[1])); else err[1]=0.; + if(fabs(f[2])>SMALL) err[2]=fabs(f[2])/(1.e-20*uu[UU]+fabs(uu[VY])+fabs(uu0[VY])+fabs(dt*gdetu*Gi[2])); else err[2]=0.; + if(fabs(f[3])>SMALL) err[3]=fabs(f[3])/(1.e-20*uu[UU]+fabs(uu[VZ])+fabs(uu0[VZ])+fabs(dt*gdetu*Gi[3])); else err[3]=0.; } else if(whichprim==RAD) //rad-primitives { @@ -2711,9 +2711,6 @@ calc_all_Gi_with_state(ldouble *pp, void *sss, void* ggg, ldouble Ruu = Ehatrad; //Calculate G thermal in lab frame -#ifdef APPLY_OMP_SIMD - //#pragma omp simd -#endif for(i = 0; i < 4; i++) { ldouble Ru = 0.; @@ -2812,9 +2809,9 @@ calc_Gi_nonrel_with_state(ldouble *pp, void *sss, void *ggg, ldouble Gi[4], int //the four-velocity of fluid in lab frame ldouble ucon[4],utcon[4],ucov[4],vpr[3],vgas[3]; - utcon[1]=pp[2]; - utcon[2]=pp[3]; - utcon[3]=pp[4]; + utcon[1]=pp[VX]; + utcon[2]=pp[VY]; + utcon[3]=pp[VZ]; conv_vels_both(utcon,ucon,ucov,VELPRIM,VEL4,gg,GG); vgas[0]=utcon[1]; vgas[1]=utcon[2]; @@ -2822,7 +2819,7 @@ calc_Gi_nonrel_with_state(ldouble *pp, void *sss, void *ggg, ldouble Gi[4], int //gas properties ldouble rho=pp[RHO]; - ldouble u=pp[1]; + ldouble u=pp[UU]; ldouble p= (gamma-1.)*(ldouble)u; ldouble Ti,Te; ldouble Tgas=calc_PEQ_Teifrompp(pp,&Te,&Ti,geom->ix,geom->iy,geom->iz); @@ -3071,9 +3068,6 @@ calc_Ehat_from_Rij_ucov(double Rij[4][4], double uffcov[4], ldouble *Ehat) *Ehat = 0.; for(i=0;i<4;i++) -#ifdef APPLY_OMP_SIMD - //#pragma omp simd -#endif for(j=0;j<4;j++) *Ehat+=Rij[i][j]*uffcov[i]*uffcov[j]; @@ -3540,9 +3534,9 @@ calc_ncompt_Thatrad_full(ldouble *pp, void* ggg) //the four-velocity of fluid ldouble uffcon[4],utcon[4],uffcov[4],vpr[3]; - utcon[1]=pp[2]; - utcon[2]=pp[3]; - utcon[3]=pp[4]; + utcon[1]=pp[VX]; + utcon[2]=pp[VY]; + utcon[3]=pp[VZ]; conv_vels_both(utcon,uffcon,uffcov,VELPRIM,VEL4,gg,GG); //R^ab u_a u_b = Erad in fluid frame @@ -3580,9 +3574,9 @@ calc_ncompt_Thatrad(ldouble *pp, void* ggg, ldouble Ehatrad) //the four-velocity of fluid in lab frame ldouble uffcon[4],utcon[4],uffcov[4],vpr[3]; - utcon[1]=pp[2]; - utcon[2]=pp[3]; - utcon[3]=pp[4]; + utcon[1]=pp[VX]; + utcon[2]=pp[VY]; + utcon[3]=pp[VZ]; conv_vels_both(utcon,uffcon,uffcov,VELPRIM,VEL4,gg,GG); //the four-velocity of radiation in lab frame @@ -3688,9 +3682,9 @@ calc_ncompt_nphhat(ldouble *pp, void* ggg) //the four-velocity of fluid in lab frame ldouble uffcon[4],utcon[4],uffcov[4],vpr[3]; - utcon[1]=pp[2]; - utcon[2]=pp[3]; - utcon[3]=pp[4]; + utcon[1]=pp[VX]; + utcon[2]=pp[VY]; + utcon[3]=pp[VZ]; conv_vels_both(utcon,uffcon,uffcov,VELPRIM,VEL4,gg,GG); //the four-velocity of radiation in lab frame @@ -4806,9 +4800,6 @@ calc_Rij_visc(ldouble *pp, void* ggg, ldouble Rvisc[][4], int *derdir) for(i=0;i<4;i++) { -#ifdef APPLY_OMP_SIMD - //#pragma omp simd -#endif for(j=0;j<4;j++) { Rvisc[i][j]=0.; @@ -5256,15 +5247,15 @@ test_solve_implicit_lab() struct geometry geom; fill_geometry(TNX-1,TNY/2,0,&geom); - pp0[0]=1.e-21; - pp0[1]=1.e-23;//calc_PEQ_ufromTrho(1.e9,pp0[0]); - pp0[2]=0.00001; - pp0[3]=0.0; - pp0[4]=0.0; - pp0[5]=calc_Sfromu(pp0[0],pp0[1],TNX-1,TNY/2,0); + pp0[RHO]=1.e-21; + pp0[UU]=1.e-23; + pp0[VX]=0.00001; + pp0[VY]=0.0; + pp0[VZ]=0.0; + pp0[ENTR]=calc_Sfromu(pp0[RHO],pp0[UU],TNX-1,TNY/2,0); pp0[B1]=pp0[B2]=pp0[B3]=0.; - pp0[EE]=1.e-23;//calc_LTE_Efromurho(pp0[1],pp0[0]);//1.e-21; + pp0[EE]=1.e-23; pp0[FX]=0.00001; pp0[FY]=0.0; pp0[FZ]=0.0; @@ -5287,7 +5278,7 @@ test_solve_implicit_lab() PLOOP(iv) pp[iv]=pp0[iv]; print_primitives(pp0); - printf("gas temp: %e\nrad temp: %e\n",calc_PEQ_Tfromurho(pp0[1],pp0[0],geom.ix,geom.iy,geom.iz),calc_LTE_TfromE(pp0[EE])); + printf("gas temp: %e\nrad temp: %e\n",calc_PEQ_Tfromurho(pp0[UU],pp0[RHO],geom.ix,geom.iy,geom.iz),calc_LTE_TfromE(pp0[EE])); ldouble del4[NRADVAR]; int verbose=0; @@ -5331,7 +5322,6 @@ test_solve_implicit_lab() PLOOP(iv) pp[iv]=pp0[iv]; if(1) { - params[0]=MHD; params[1]=RADIMPLICIT_ENERGYEQ; params[2]=RADIMPLICIT_LAB; @@ -5340,7 +5330,7 @@ test_solve_implicit_lab() params[6]=0; //ifelectron solve_implicit_lab_4dprim(uu0,pp0,&geom,dt,verbose,params,pp); print_primitives(pp); - //printf("gas temp: %e\nrad temp: %e\n",calc_PEQ_Tfromurho(pp[1],pp[0]),calc_LTE_TfromE(pp[EE])); + printf("gas temp: %e\nrad temp: %e\n",calc_PEQ_Tfromurho(pp[UU],pp[RHO]),calc_LTE_TfromE(pp[EE])); } #endif exit(1); @@ -5379,12 +5369,12 @@ test_Gi() printf("Giff[0] = %.3e\n",Gic[0]*conv); exit(1); - pp0[0]=9.253e-9; - pp0[1]=calc_PEQ_ufromTrho(2.478e9,pp0[RHO],geom.ix,geom.iy,geom.iz); + pp0[RHO]=9.253e-9; + pp0[UU]=calc_PEQ_ufromTrho(2.478e9,pp0[RHO],geom.ix,geom.iy,geom.iz); pp0[VX]=-1.04792e-01; pp0[VY]=2.69291e-03; pp0[VZ]=1.62996e-02; - pp0[5]=calc_Sfromu(pp0[0],pp0[1],geom.ix,geom.iy,geom.iz); + pp0[ENTR]=calc_Sfromu(pp0[RHO],pp0[UU],geom.ix,geom.iy,geom.iz); pp0[EE]=1.0*calc_LTE_Efromurho(pp0[UU],pp0[RHO]); pp0[FX]=0.01; pp0[FY]=0.; @@ -5457,9 +5447,9 @@ test_solve_implicit_lab_file() ldouble vcon[4],ucon[4],ucov[4]; //converting to 4-velocity - vcon[1]=pp0[2]; - vcon[2]=pp0[3]; - vcon[3]=pp0[4]; + vcon[1]=pp0[VX]; + vcon[2]=pp0[VY]; + vcon[3]=pp0[VZ]; vcon[0]=0.; conv_vels_both(vcon,ucon,ucov,VELPRIM,VEL4,geom.gg,geom.GG); print_4vector(ucon); @@ -5533,7 +5523,7 @@ test_jon_solve_implicit_lab() ucon[3]=pp[VZ]; conv_vels(ucon,ucon,VEL4,VEL4,geom.gg,geom.GG); geom.alpha=sqrt(-1./geom.GG[0][0]); - pp[5]=calc_Sfromu(pp[0],pp[1],geom.ix,geom.iy,geom.iz); + pp[ENTR]=calc_Sfromu(pp[RHO],pp[UU],geom.ix,geom.iy,geom.iz); //destroy magn field //uu[B1]=uu[B2]=uu[B3]=pp[B1]=pp[B2]=pp[B3]=0.; diff --git a/relele.c b/relele.c index 1d83558..5fed3ed 100644 --- a/relele.c +++ b/relele.c @@ -313,9 +313,6 @@ calc_alpgam(ldouble *u1, ldouble gg[][5], ldouble GG[][5]) int i, j; ldouble qsq=0.; -#ifdef APPLY_OMP_SIMD - //#pragma omp simd -#endif for(i=1;i<4;i++) { for(j=1;j<4;j++) @@ -422,15 +419,15 @@ conv_velsinprims(ldouble *pp,int which1, int which2,ldouble gg[][5],ldouble GG[] int ret=0; ldouble v1[4],v2[4]; v1[0]=0.; //not considered - v1[1]=pp[2]; - v1[2]=pp[3]; - v1[3]=pp[4]; + v1[1]=pp[VX]; + v1[2]=pp[VY]; + v1[3]=pp[VZ]; ret=conv_vels(v1,v2,which1,which2,gg,GG); if(ret==0) { - pp[2]=v2[1]; - pp[3]=v2[2]; - pp[4]=v2[3]; + pp[VX]=v2[1]; + pp[VY]=v2[2]; + pp[VZ]=v2[3]; return 0; } else @@ -539,16 +536,16 @@ set_hdatmosphere(ldouble *pp,ldouble *xx,ldouble gg[][5],ldouble GG[][5],int atm // to VELPRIM conv_vels(ucon,ucon,VEL4,VELPRIM,gg,GG); - pp[2]=ucon[1]; - pp[3]=ucon[2]; - pp[4]=ucon[3]; + pp[VX]=ucon[1]; + pp[VY]=ucon[2]; + pp[VZ]=ucon[3]; // density & pressure ldouble r=xx2[1]; ldouble rout=1.; //RHOATMMIN etc. given at rout=2 - pp[0] = RHOATMMIN*pow(r/rout,-1.5); - pp[1] = UINTATMMIN*pow(r/rout,-2.5); + pp[RHO] = RHOATMMIN*pow(r/rout,-1.5); + pp[UU] = UINTATMMIN*pow(r/rout,-2.5); #ifdef MAGNFIELD pp[B1]=pp[B2]=pp[B3]=0.; @@ -566,9 +563,9 @@ set_hdatmosphere(ldouble *pp,ldouble *xx,ldouble gg[][5],ldouble GG[][5],int atm { conv_vels(ucon,ucon,VELR,VELPRIM,gg,GG); } - pp[2]=ucon[1]; - pp[3]=ucon[2]; - pp[4]=ucon[3]; + pp[VX]=ucon[1]; + pp[VY]=ucon[2]; + pp[VZ]=ucon[3]; // Bondi-like atmosphere coco_N(xx,xx2,MYCOORDS,BLCOORDS); @@ -576,8 +573,8 @@ set_hdatmosphere(ldouble *pp,ldouble *xx,ldouble gg[][5],ldouble GG[][5],int atm ldouble rout=2.; //RHOATMMIN etc. given at rout=2 - pp[0] = RHOATMMIN*pow(r/rout,-1.5); - pp[1] = UINTATMMIN*pow(r/rout,-2.5); + pp[RHO] = RHOATMMIN*pow(r/rout,-1.5); + pp[UU] = UINTATMMIN*pow(r/rout,-2.5); #ifdef MAGNFIELD pp[B1]=pp[B2]=pp[B3]=0.; @@ -595,9 +592,9 @@ set_hdatmosphere(ldouble *pp,ldouble *xx,ldouble gg[][5],ldouble GG[][5],int atm { conv_vels(ucon,ucon,VELR,VELPRIM,gg,GG); } - pp[2]=ucon[1]; - pp[3]=ucon[2]; - pp[4]=ucon[3]; + pp[VX]=ucon[1]; + pp[VY]=ucon[2]; + pp[VZ]=ucon[3]; // Bondi-like atmosphere coco_N(xx,xx2,MYCOORDS,BLCOORDS); @@ -605,8 +602,8 @@ set_hdatmosphere(ldouble *pp,ldouble *xx,ldouble gg[][5],ldouble GG[][5],int atm ldouble rout=2.; //RHOATMMIN etc. given at rout=2 - pp[0] = RHOATMMIN*pow(r/rout,-2.0); - pp[1] = UINTATMMIN*pow(r/rout,-2.5); + pp[RHO] = RHOATMMIN*pow(r/rout,-2.0); + pp[UU] = UINTATMMIN*pow(r/rout,-2.5); #ifdef MAGNFIELD pp[B1]=pp[B2]=pp[B3]=0.; @@ -624,11 +621,11 @@ set_hdatmosphere(ldouble *pp,ldouble *xx,ldouble gg[][5],ldouble GG[][5],int atm { conv_vels(ucon,ucon,VELR,VELPRIM,gg,GG); } - pp[2]=ucon[1]; - pp[3]=ucon[2]; - pp[4]=ucon[3]; - pp[0] = RHOATMMIN; - pp[1] = UINTATMMIN; + pp[VX]=ucon[1]; + pp[VY]=ucon[2]; + pp[VZ]=ucon[3]; + pp[RHO] = RHOATMMIN; + pp[UU] = UINTATMMIN; #ifdef MAGNFIELD pp[B1]=pp[B2]=pp[B3]=0.; @@ -657,7 +654,10 @@ set_hdatmosphere(ldouble *pp,ldouble *xx,ldouble gg[][5],ldouble GG[][5],int atm //corrected rho: rho=PAR_D/(r*r*sqrtl(2./r)); - pp[0]=rho; pp[1]=uint; pp[2]=-V; pp[3]=pp[4]=0.; + pp[RHO]=rho; + pp[UU]=uint; + pp[VX]=-V; + pp[VY]=pp[VZ]=0.; conv_velsinprims(pp,VEL3,VELPRIM,gg,GG); #ifdef MAGNFIELD @@ -686,17 +686,17 @@ set_hdatmosphere(ldouble *pp,ldouble *xx,ldouble gg[][5],ldouble GG[][5],int atm // to VELPRIM conv_vels(ucon,ucon,VEL4,VELPRIM,gg,GG); - pp[2]=ucon[1]; - pp[3]=ucon[2]; - pp[4]=ucon[3]; + pp[VX]=ucon[1]; + pp[VY]=ucon[2]; + pp[VZ]=ucon[3]; //density etc. ldouble r=xx2[1]; ldouble rout=2.; //RHOATMMIN etc. given at rout=2 - pp[0] = RHOATMMIN*pow(r/rout,-1.5); - pp[1] = UINTATMMIN*pow(r/rout,-2.5); + pp[RHO] = RHOATMMIN*pow(r/rout,-1.5); + pp[UU] = UINTATMMIN*pow(r/rout,-2.5); #ifdef MAGNFIELD pp[B1]=pp[B2]=pp[B3]=0.; @@ -1057,7 +1057,7 @@ pick_Gb(int ix,int iy,int iz,int idim,ldouble gg[][5]) int print_p(ldouble *p) { - printf("rho: %10e\nuu: %10e\nvr: %10e\nvth: %10e\nvph: %10e\nS: %10e\n",p[0],p[1],p[2],p[3],p[4],p[5]); + printf("rho: %10e\nuu: %10e\nvr: %10e\nvth: %10e\nvph: %10e\nS: %10e\n",p[RHO],p[UU],p[VX],p[VY],p[VZ],p[ENTR]); #ifdef RADIATION printf("E: %10e\nFx: %10e\nFy: %10e\nFz: %10e\n\n",p[EE0],p[FX0],p[FY0],p[FZ0]); #endif @@ -1072,7 +1072,7 @@ print_p(ldouble *p) int print_u(ldouble *u) { - printf("rhout: %10e\nTtt: %10e\nTtr: %10e\nTtth: %10e\nTtph: %10e\nSut: %10e\n",u[0],u[1]-u[0],u[2],u[3],u[4],u[5]); + printf("rhout: %10e\nTtt: %10e\nTtr: %10e\nTtth: %10e\nTtph: %10e\nSut: %10e\n",u[RHO],u[UU]-u[RHO],u[VX],u[VY],u[VZ],u[ENTR]); #ifdef RADIATION printf("Rtt: %10e\nRt1: %10e\nRt2: %10e\nRt3: %10e\n\n",u[EE0],u[FX0],u[FY0],u[FZ0]); #endif diff --git a/silo.c b/silo.c index e77763e..6760b7f 100644 --- a/silo.c +++ b/silo.c @@ -24,9 +24,9 @@ int fprint_silofile(ldouble time, int num, char* folder, char* prefix) mpi_exchangedata(); calc_avgs_throughout(); - DBfile *file = NULL; // The Silo file pointer - char *coordnames[3]; // Names of the coordinates - ldouble *nodex; // The coordinate arrays + DBfile *file = NULL; // The Silo file pointer + char *coordnames[3]; // Names of the coordinates + ldouble *nodex; // The coordinate arrays ldouble *nodey; ldouble *nodez; ldouble *coordinates[3];// The array of coordinatearrays @@ -60,13 +60,6 @@ int fprint_silofile(ldouble time, int num, char* folder, char* prefix) nz++;nnz++; #endif - //initialize shared gammabreak array for nonthermal electrons -#ifdef RELELECTRONS - int ie; - ldouble gammapbrk[NRELBIN]; - for(ie=0; ie0) continue; - if(NY==1 && iy>0) continue; - - int iix,iiy,iiz; - iix=ix; - iiy=iy; - iiz=iz; - - if(NZ>1 && NY>1) - { - xxvec[0]=0.;xxvec[1]=get_xb(iix,0);xxvec[2]=get_xb(iiy,1);xxvec[3]=get_xb(iiz,2); - } - else if(NZ==1 && NY==1) - { - xxvec[0]=0.;xxvec[1]=get_xb(iix,0);xxvec[2]=get_x(iiy,1);xxvec[3]=get_x(iiz,2); - } - else if(NZ==1) - { - xxvec[0]=0.;xxvec[1]=get_xb(iix,0);xxvec[2]=get_xb(iiy,1);xxvec[3]=get_x(iiz,2); - } - else if(NY==1) - { - xxvec[0]=0.;xxvec[1]=get_xb(iix,0);xxvec[2]=get_x(iiy,1);xxvec[3]=get_xb(iiz,2); - } - - coco_N(xxvec,xxvecsph,MYCOORDS,SPHCOORDS); - coco_N(xxvec,xxveccar,MYCOORDS,MINKCOORDS); - - imz=iz; - imy=iy; - imx=ix; -#ifdef PRINTZGC_RIGHT - imz=iz-NG; -#endif -#ifdef PRINTXGC_RIGHT - imx=ix-NG; -#endif -#ifdef PRINTXGC_LEFT - imx=ix+NG; -#endif -#ifdef PRINTYGC_RIGHT - imy=iy-NG; -#endif + if(NZ==1 && iz>0) continue; + if(NY==1 && iy>0) continue; + + int iix,iiy,iiz; + iix=ix; + iiy=iy; + iiz=iz; + + if(NZ>1 && NY>1) + { + xxvec[0]=0.;xxvec[1]=get_xb(iix,0);xxvec[2]=get_xb(iiy,1);xxvec[3]=get_xb(iiz,2); + } + else if(NZ==1 && NY==1) + { + xxvec[0]=0.;xxvec[1]=get_xb(iix,0);xxvec[2]=get_x(iiy,1);xxvec[3]=get_x(iiz,2); + } + else if(NZ==1) + { + xxvec[0]=0.;xxvec[1]=get_xb(iix,0);xxvec[2]=get_xb(iiy,1);xxvec[3]=get_x(iiz,2); + } + else if(NY==1) + { + xxvec[0]=0.;xxvec[1]=get_xb(iix,0);xxvec[2]=get_x(iiy,1);xxvec[3]=get_xb(iiz,2); + } + + coco_N(xxvec,xxvecsph,MYCOORDS,SPHCOORDS); + coco_N(xxvec,xxveccar,MYCOORDS,MINKCOORDS); + + imz=iz; + imy=iy; + imx=ix; + #ifdef PRINTZGC_RIGHT + imz=iz-NG; + #endif + #ifdef PRINTXGC_RIGHT + imx=ix-NG; + #endif + #ifdef PRINTXGC_LEFT + imx=ix+NG; + #endif + #ifdef PRINTYGC_RIGHT + imy=iy-NG; + #endif + + int nodalindex=imz*(nny*nnx) + imy*nnx + imx; + + //coordinates + + ldouble coordscale=1.; + #ifdef COORDSINPCINSILO + coordscale=1./(PARSECCGS/MASSCM); + #endif - int nodalindex=imz*(nny*nnx) + imy*nnx + imx; - - //coordinates - - ldouble coordscale=1.; - #ifdef COORDSINPCINSILO - coordscale=1./(PARSECCGS/MASSCM); - #endif - #if(SILOCOORDS==SPHCOORDS) - nodex[nodalindex]=xxvecsph[1]*coordscale; - nodey[nodalindex]=xxvecsph[2]*coordscale; - nodez[nodalindex]=xxvecsph[3]*coordscale; + nodex[nodalindex]=xxvecsph[1]*coordscale; + nodey[nodalindex]=xxvecsph[2]*coordscale; + nodez[nodalindex]=xxvecsph[3]*coordscale; #else - nodex[nodalindex]=xxveccar[1]*coordscale; - nodey[nodalindex]=xxveccar[2]*coordscale; - nodez[nodalindex]=xxveccar[3]*coordscale; -#endif - - - - } - } -} + nodex[nodalindex]=xxveccar[1]*coordscale; + nodey[nodalindex]=xxveccar[2]*coordscale; + nodez[nodalindex]=xxveccar[3]*coordscale; +#endif + } + } + } //then fill the zones with values @@ -302,954 +304,910 @@ int fprint_silofile(ldouble time, int num, char* folder, char* prefix) #else for(iz=0;iz1 && iy==0) || (NZ>1 && iz==0)) //divB left-biased - divB[zonalindex]=0; - else - divB[zonalindex]=calc_divB(ix,iy,iz); - #endif - - taucoupling[zonalindex]=estimate_gas_radiation_coupling(pp,&geomout); - - - dpdr = (gdet2*gamma*get_u(p,UU,iix+1,iiy,iiz)-gdet1*gamma*get_u(p,UU,iix-1,iiy,iiz)) / (xx2[1]-xx1[1]); - gracen=0.; - for(i=0;i<4;i++) - for(j=0;j<4;j++) - gracen += gdet*Tij[i][j]*get_gKr(j,1,i,ix,iy,iz); - - - #ifdef PRINTVISCHEATINGTOSILO - deltae[zonalindex]=calc_ViscousElectronHeatingFraction(&get_u(p,0,ix,iy,iz),&geomout); - vischeat[zonalindex]=get_u_scalar(vischeating,ix,iy,iz); - vischeatnege[zonalindex]=get_u_scalar(vischeatingnegebalance,ix,iy,iz);; - vischeatnegi[zonalindex]=get_u_scalar(vischeatingnegibalance,ix,iy,iz);; - dtauarr[zonalindex]=-1.; - #endif - } - - else //using averaged data - { - entr[zonalindex]=get_uavg(pavg,ENTR,ix,iy,iz); - rho[zonalindex]=get_uavg(pavg,RHO,ix,iy,iz); - uint[zonalindex]=get_uavg(pavg,UU,ix,iy,iz); - - pregas=get_uavg(pavg,AVGPGAS,ix,iy,iz); - gamma = 1. + pregas/uint[zonalindex]; - - #ifdef RADIATION - prerad=(1./3.)*get_uavg(pavg,AVGEHAT,ix,iy,iz); - #else - prerad=0.; - #endif - - #ifdef MAGNFIELD - premag=0.5*bsq[zonalindex]; - //beta[zonalindex]=(pregas + prerad + premag)/premag; - beta[zonalindex]=pregas/0.5/bsq[zonalindex]; - betainv[zonalindex]=1./beta[zonalindex]; - sigma[zonalindex]=bsq[zonalindex]/(rho[zonalindex] + uint[zonalindex] + pregas); - #endif - - rhouconr=get_uavg(pavg,AVGRHOUCON(1),ix,iy,iz); - rhoucont=get_uavg(pavg,AVGRHOUCON(0),ix,iy,iz); - - vel[0]=get_uavg(pavg,AVGRHOUCON(0),ix,iy,iz)/get_uavg(pavg,RHO,ix,iy,iz); - vel[1]=get_uavg(pavg,AVGRHOUCON(1),ix,iy,iz)/get_uavg(pavg,RHO,ix,iy,iz); - vel[2]=get_uavg(pavg,AVGRHOUCON(2),ix,iy,iz)/get_uavg(pavg,RHO,ix,iy,iz); - vel[3]=get_uavg(pavg,AVGRHOUCON(3),ix,iy,iz)/get_uavg(pavg,RHO,ix,iy,iz); - for(i=0;i<4;i++) vcon[i]=vel[i]; - indices_21(vel,vcov,geomout.gg); - - Omega[zonalindex]=vel[3]/vel[0]; - pp[VX]=vel[1]; //updates pp[VI] to have rho-weighted velocities there - pp[VY]=vel[2]; - pp[VZ]=vel[3]; - - for(i=0;i<4;i++) - for(j=0;j<4;j++) - Tij[i][j]=get_uavg(pavg,AVGRHOUCONUCOV(i,j),ix,iy,iz) - + GAMMA*get_uavg(pavg,AVGUUUCONUCOV(i,j),ix,iy,iz) - + get_uavg(pavg,AVGBSQUCONUCOV(i,j),ix,iy,iz) - + delta(i,j)*(GAMMA*get_uavg(pavg,UU,ix,iy,iz) + 1./2.*get_uavg(pavg,AVGBSQ,ix,iy,iz)) - - get_uavg(pavg,AVGBCONBCOV(i,j),ix,iy,iz); - - //overwrite with directly averaged Tij - for(i=0;i<4;i++) - for(j=0;j<4;j++) - Tij[i][j]=get_uavg(pavg,AVGTIJ(i,j),ix,iy,iz); - - - Tit[1]=Tij[1][0]; - Tit[2]=Tij[2][0]; - Tit[3]=Tij[3][0]; - - //Bernoulli flux / number - muBe[zonalindex]=-(Tij[1][0] + rhouconr)/rhouconr; - Be[zonalindex]=-(Tij[0][0] + rhoucont)/rhoucont; - #ifdef MAGNFIELD - - Qtheta[zonalindex]=2.*M_PI/Omega[zonalindex]/dx[1]*fabs(bcon[2])/sqrt(rho[zonalindex]); - Qphi[zonalindex]=2.*M_PI/Omega[zonalindex]/dx[2]*fabs(bcon[3])/sqrt(rho[zonalindex]); - ldouble brbphi,bsq,bfake[4]; - - #ifdef BHDISK_PROBLEMTYPE - calc_angle_brbphibsq(ix,iy,iz,&brbphi,&bsq,bfake,bfake); //to calculate magn. field angle - Bangle[zonalindex]=-brbphi/bsq; - #else - Bangle[zonalindex]=-1.; - #endif - - if(ix==0 || (NY>1 && iy==0) || (NZ>1 && iz==0)) //divB left-biased - divB[zonalindex]=0; - else - divB[zonalindex]=calc_divB(ix,iy,iz); - #endif - - taucoupling[zonalindex]=estimate_gas_radiation_coupling(pp,&geomout); - - dpdr = (gdet2*GAMMA*get_uavg(pavg,UU,iix+1,iiy,iiz)-gdet1*GAMMA*get_uavg(pavg,UU,iix-1,iiy,iiz)) / (xx2[1]-xx1[1]); - gracen=0.; - for(i=0;i<4;i++) - for(j=0;j<4;j++) - gracen += gdet*Tij[i][j]*get_gKr(j,1,i,ix,iy,iz); - - #ifdef EVOLVEELECTRONS - ldouble pe,pi; - pe=get_uavg(pavg,AVGPE,ix,iy,iz); - pi=get_uavg(pavg,AVGPI,ix,iy,iz); - #ifndef CONSISTENTGAMMA - pi=pregas-pe; - #endif - - //electrons - ne=calc_thermal_ne(pp); - tempeloc=pe/K_BOLTZ/ne; - ldouble gammae=GAMMAE; - #ifdef CONSISTENTGAMMA - #ifndef FIXEDGAMMASPECIES - gammae=calc_gammaintfromtemp(tempeloc,ELECTRONS); - #endif - #endif - ueloc=pe/(gammae-1.); - - //ions - ldouble ni=rho[zonalindex]/MU_I/M_PROTON; - tempiloc=pi/K_BOLTZ/ni; - ldouble gammai=GAMMAI; - #ifdef CONSISTENTGAMMA - #ifndef FIXEDGAMMASPECIES - gammai=calc_gammaintfromtemp(tempiloc,IONS); - #endif - #endif - uiloc=pi/(gammai-1.); - #endif // EVOLVEELECTRONS + gdet1=calc_gdet_arb(xx1,OUTCOORDS); + gdet2=calc_gdet_arb(xx2,OUTCOORDS); + + // compute zonal index + imz=iz; + imy=iy; + imx=ix; + #ifdef PRINTZGC_RIGHT + imz=iz-NG; + #endif + #ifdef PRINTXGC_RIGHT + imx=ix-NG; + #endif + #ifdef PRINTXGC_LEFT + imx=ix+NG; + #endif + #ifdef PRINTYGC_RIGHT + imy=iy-NG; + #endif + + int zonalindex=imz*(ny*nx) + imy*nx + imx; + + // copy over primitives + for(iv=0;iv1 && iy==0) || (NZ>1 && iz==0)) //divB left-biased + divB[zonalindex]=0; + else + divB[zonalindex]=calc_divB(ix,iy,iz); +#endif //MAGNFIELD + +#ifdef RADIATION + taucoupling[zonalindex]=estimate_gas_radiation_coupling(pp,&geomout); +#endif //RADIATION + + // electrons & ions + ueloc=uiloc=0.; +#ifdef EVOLVEELECTRONS + // temperatures + calc_PEQ_Teifrompp(pp,&tempeloc,&tempiloc,ix,iy,iz); + + // electrons + ldouble ne=calc_thermal_ne(pp); + ldouble pe=K_BOLTZ*ne*tempeloc; + ldouble gammae=GAMMAE; + #ifdef CONSISTENTGAMMA + #ifndef FIXEDGAMMASPECIES + gammae=calc_gammaintfromtemp(tempeloc,ELECTRONS); + #endif + #endif + ueloc=pe/(gammae-1.); + + //ions + ldouble ni=pp[RHO]/MU_I/M_PROTON; //number density of photons and electrons + ldouble pi=K_BOLTZ*ni*tempiloc; + ldouble gammai=GAMMAI; + #ifdef CONSISTENTGAMMA + #ifndef FIXEDGAMMASPECIES + gammai=calc_gammaintfromtemp(tempiloc,IONS); + #endif + #endif + uiloc=pi/(gammai-1.); +#endif //EVOLVEELECTRONS #ifdef PRINTVISCHEATINGTOSILO - vischeat[zonalindex]=get_uavg(pavg,AVGVISCHEATING,ix,iy,iz); - vischeatnege[zonalindex]=get_uavg(pavg,AVGVISCHEATINGNEGE,ix,iy,iz); - vischeatnegi[zonalindex]=get_uavg(pavg,AVGVISCHEATINGNEGI,ix,iy,iz); - deltae[zonalindex]=-1.; - dtauarr[zonalindex]=-1.; -#endif - } //doingavg + deltae[zonalindex]=calc_ViscousElectronHeatingFraction(&get_u(p,0,ix,iy,iz),&geomout); + vischeat[zonalindex]=get_u_scalar(vischeating,ix,iy,iz); + vischeatnege[zonalindex]=get_u_scalar(vischeatingnegebalance,ix,iy,iz);; + vischeatnegi[zonalindex]=get_u_scalar(vischeatingnegibalance,ix,iy,iz);; + dtauarr[zonalindex]=-1.; +#endif //PRINTVISCHEATINGTOSILO + } + else //using averaged data + { + rho[zonalindex] = get_uavg(pavg,RHO,ix,iy,iz); + entr[zonalindex] = get_uavg(pavg,ENTR,ix,iy,iz); + uint[zonalindex] = get_uavg(pavg,UU,ix,iy,iz); + + pregas = get_uavg(pavg,AVGPGAS,ix,iy,iz); + gamma = 1. + pregas/uint[zonalindex]; + + // velocity + vel[0]=get_uavg(pavg,AVGRHOUCON(0),ix,iy,iz)/get_uavg(pavg,RHO,ix,iy,iz); + vel[1]=get_uavg(pavg,AVGRHOUCON(1),ix,iy,iz)/get_uavg(pavg,RHO,ix,iy,iz); + vel[2]=get_uavg(pavg,AVGRHOUCON(2),ix,iy,iz)/get_uavg(pavg,RHO,ix,iy,iz); + vel[3]=get_uavg(pavg,AVGRHOUCON(3),ix,iy,iz)/get_uavg(pavg,RHO,ix,iy,iz); + for(i=0;i<4;i++) vcon[i]=vel[i]; + indices_21(vcon,vcov,geomout.gg); + + rhouconr=get_uavg(pavg,AVGRHOUCON(1),ix,iy,iz); + rhoucont=get_uavg(pavg,AVGRHOUCON(0),ix,iy,iz); + + // angular velocity + Omega[zonalindex]=vel[3]/vel[0]; + + // ANDREW these are 4-velocites, not VELR, so I don't think this is right... + //updates pp[VI] to have rho-weighted velocities there + //pp[VX]=vel[1]; + //pp[VY]=vel[2]; + //pp[VZ]=vel[3]; + + #ifdef MAGNFIELD + B1comp[zonalindex] = get_uavg(pavg,B1,ix,iy,iz); + B2comp[zonalindex] = get_uavg(pavg,B2,ix,iy,iz); + B3comp[zonalindex] = get_uavg(pavg,B3,ix,iy,iz); + + bcon[1]=get_uavg(pavg,AVGBCON(1),ix,iy,iz); + bcon[2]=get_uavg(pavg,AVGBCON(2),ix,iy,iz); + bcon[3]=get_uavg(pavg,AVGBCON(3),ix,iy,iz); + bsq[zonalindex]=get_uavg(pavg,AVGBSQ,ix,iy,iz); + + premag=0.5*bsq[zonalindex]; + beta[zonalindex]=pregas/0.5/bsq[zonalindex]; + betainv[zonalindex]=1./beta[zonalindex]; + sigma[zonalindex]=bsq[zonalindex]/(rho[zonalindex] + uint[zonalindex] + pregas); + sigmaw[zonalindex]=bsq[zonalindex]/(rho[zonalindex] + gamma*uint[zonalindex]); + + // field line angular velocity + //OmegaF[zonalindex] = vel[3]/vel[0] - (vel[1]/vel[0])*(pp[B3]/pp[B1]); + OmegaF[zonalindex] = (Omega[zonalindex] - (vel[2]*B3comp[zonalindex])/(vel[0]*B2comp[zonalindex])); + #endif + + // stress energy tensor + /* + for(i=0;i<4;i++) + for(j=0;j<4;j++) + Tij[i][j]=get_uavg(pavg,AVGRHOUCONUCOV(i,j),ix,iy,iz) + + GAMMA*get_uavg(pavg,AVGUUUCONUCOV(i,j),ix,iy,iz) + + get_uavg(pavg,AVGBSQUCONUCOV(i,j),ix,iy,iz) + + delta(i,j)*(GAMMA*get_uavg(pavg,UU,ix,iy,iz) + 1./2.*get_uavg(pavg,AVGBSQ,ix,iy,iz)) + - get_uavg(pavg,AVGBCONBCOV(i,j),ix,iy,iz); + */ + + // averaged stress energy tensor + for(i=0;i<4;i++) + for(j=0;j<4;j++) + Tij[i][j]=get_uavg(pavg,AVGTIJ(i,j),ix,iy,iz); + + + Tit[1]=Tij[1][0]; + Tit[2]=Tij[2][0]; + Tit[3]=Tij[3][0]; + + // Bernoulli flux / number + muBe[zonalindex]=-(Tij[1][0] + rhouconr)/rhouconr; + Be[zonalindex]=-(Tij[0][0] + rhoucont)/rhoucont; + +#ifdef MAGNFIELD + + Qtheta[zonalindex]=2.*M_PI/Omega[zonalindex]/dx[1]*fabs(bcon[2])/sqrt(rho[zonalindex]); + Qphi[zonalindex]=2.*M_PI/Omega[zonalindex]/dx[2]*fabs(bcon[3])/sqrt(rho[zonalindex]); + ldouble brbphi,bsq,bfake[4]; + + #ifdef BHDISK_PROBLEMTYPE + calc_angle_brbphibsq(ix,iy,iz,&brbphi,&bsq,bfake,bfake); //to calculate magn. field angle + Bangle[zonalindex]=-brbphi/bsq; + #else + Bangle[zonalindex]=-1.; + #endif + + if(ix==0 || (NY>1 && iy==0) || (NZ>1 && iz==0)) //divB left-biased + divB[zonalindex]=0; + else + divB[zonalindex]=calc_divB(ix,iy,iz); +#endif // MAGNFIELD + +#ifdef RADIATION + taucoupling[zonalindex]=estimate_gas_radiation_coupling(pp,&geomout); +#endif // RADIATION + +#ifdef EVOLVEELECTRONS + ldouble pe,pi; + pe=get_uavg(pavg,AVGPE,ix,iy,iz); + pi=get_uavg(pavg,AVGPI,ix,iy,iz); + #ifndef CONSISTENTGAMMA + pi=pregas-pe; + #endif + + //electrons + ne=calc_thermal_ne(pp); + tempeloc=pe/K_BOLTZ/ne; + ldouble gammae=GAMMAE; + #ifdef CONSISTENTGAMMA + #ifndef FIXEDGAMMASPECIES + gammae=calc_gammaintfromtemp(tempeloc,ELECTRONS); + #endif + #endif + ueloc=pe/(gammae-1.); + + //ions + ldouble ni=rho[zonalindex]/MU_I/M_PROTON; + tempiloc=pi/K_BOLTZ/ni; + ldouble gammai=GAMMAI; + #ifdef CONSISTENTGAMMA + #ifndef FIXEDGAMMASPECIES + gammai=calc_gammaintfromtemp(tempiloc,IONS); + #endif + #endif + uiloc=pi/(gammai-1.); +#endif // EVOLVEELECTRONS - /* - //ANDREW - //in avg vischeat was averaged as du, not du/dtau - //here, recompute dt and use that as an estimate #ifdef PRINTVISCHEATINGTOSILO -#ifdef DIVIDEVISCHEATBYDT - //get new prims in original coords - ldouble xxvecout[4]; - ldouble velcoord[4]; - - coco_N(xxvec,xxvecout,MYCOORDS,OUTCOORDS); - fill_utinucon(vel,geomout.gg,geomout.GG); - trans2_coco(xxvecout,vel,velcoord,OUTCOORDS,MYCOORDS); //vel should already be 4-velocity - - dt=get_u_scalar(cell_dt,ix,iy,iz); //individual time step //?? but need min? - dt=1./tstepdenmax; - - ldouble dtau = dt/velcoord[0]; - dtauarr[zonalindex] = dtau; - vischeat[zonalindex] /=dtau; - vischeatnege[zonalindex]/=dtau; - vischeatnegi[zonalindex]/=dtau; -#endif -#endif - */ - //Nonthermal electron quantities + vischeat[zonalindex]=get_uavg(pavg,AVGVISCHEATING,ix,iy,iz); + vischeatnege[zonalindex]=get_uavg(pavg,AVGVISCHEATINGNEGE,ix,iy,iz); + vischeatnegi[zonalindex]=get_uavg(pavg,AVGVISCHEATINGNEGI,ix,iy,iz); + deltae[zonalindex]=-1.; + dtauarr[zonalindex]=-1.; +#endif //PRINTVISCHEATINGTOSILO + } //doingavg + + +// Nonthermal electron quantities #ifdef EVOLVEELECTRONS #ifdef RELELECTRONS - nethloc=ne; - nrelelloc = calc_relel_ne(pp); - urelelloc = calc_relel_uint(pp); - G0relelloc = -1.*calc_relel_G0_fluidframe(pp,&geomout, 0.0, 0); //ANDREW - fluid frame - G0ic_relel_loc = -1*calc_relel_G0_fluidframe_direct(pp, &geomout, 3); - G0syn_relel_loc = -1*calc_relel_G0_fluidframe_direct(pp, &geomout, 1); - - //calculate synchrotron break frequency - //ANDREW speed this up by saving an array - gbrkloc=RELEL_INJ_MIN; - //absolute maximum of g^4*n for g > RELGAMMAMIN - ldouble nbrk=pp[NEREL(0)]*gammapbrk[0]; - ldouble nbrk2; - for(ie=1;ie nbrk) - { - nbrk=nbrk2; - gbrkloc=relel_gammas[ie]; - } - } - } + nethloc=ne; + nrelelloc = calc_relel_ne(pp); + urelelloc = calc_relel_uint(pp); + G0relelloc = -1.*calc_relel_G0_fluidframe(pp,&geomout, 0.0, 0); + G0ic_relel_loc = -1*calc_relel_G0_fluidframe_direct(pp, &geomout, 3); + G0syn_relel_loc = -1*calc_relel_G0_fluidframe_direct(pp, &geomout, 1); + + //calculate synchrotron break frequency + //ANDREW speed this up by saving an array + gbrkloc=RELEL_INJ_MIN; + //absolute maximum of g^4*n for g > RELGAMMAMIN + ldouble nbrk=pp[NEREL(0)]*gammapbrk[0]; + ldouble nbrk2; + for(ie=1;ie nbrk) + { + nbrk=nbrk2; + gbrkloc=relel_gammas[ie]; + } + } + } #endif //RELELECTRONS #endif //EVOLVEELECTRONS - - expansion[zonalindex]=exploc; - - #ifdef EVOLVEELECTRONS - tempi[zonalindex]=tempiloc; //ion temperature - tempe[zonalindex]=tempeloc; //electron temperature - ui[zonalindex]=uiloc; //ion energy density - ue[zonalindex]=ueloc; //electron energy density - - #ifdef RELELECTRONS - gammabrk[zonalindex]=gbrkloc; - urelel[zonalindex]=urelelloc; - nrelel[zonalindex]=nrelelloc; - neth[zonalindex]=nethloc; - uratio_tot[zonalindex]=urelelloc/uint[zonalindex]; - uratio_th[zonalindex]=urelelloc/ueloc; - G0relel[zonalindex]=G0relelloc; - G0icrelel[zonalindex]=G0ic_relel_loc; - G0synrelel[zonalindex]=G0syn_relel_loc; - #ifdef CGSOUTPUT - urelel[zonalindex]=endenGU2CGS(urelel[zonalindex]); - nrelel[zonalindex]=numdensGU2CGS(nrelel[zonalindex]); - neth[zonalindex]=numdensGU2CGS(neth[zonalindex]); - G0relel[zonalindex]=endenGU2CGS(G0relel[zonalindex])*timeCGS2GU(1.); - G0icrelel[zonalindex]=endenGU2CGS(G0icrelel[zonalindex])*timeCGS2GU(1.); - G0synrelel[zonalindex]=endenGU2CGS(G0synrelel[zonalindex])*timeCGS2GU(1.); - #endif - #endif - - #ifdef CGSOUTPUT - ui[zonalindex]=endenGU2CGS(ui[zonalindex]); - ue[zonalindex]=endenGU2CGS(ue[zonalindex]); - #endif - #endif - - ldouble temploc=calc_PEQ_Tfromurho(uint[zonalindex],rho[zonalindex],ix,iy,iz); - temp[zonalindex]=temploc; - entrlnT[zonalindex]=kB_over_mugas_mp*entr[zonalindex]/pp[RHO]; + + // put quantities into arrays + expansion[zonalindex]=exploc; + + lorentz[zonalindex]=fabs(vel[0])/sqrt(fabs(geomout.GG[0][0])); + u0[zonalindex]=fabs(vel[0]); + u1[zonalindex]=vel[1]; + u2[zonalindex]=vel[2]; + u3[zonalindex]=vel[3]; + + #ifdef EVOLVEELECTRONS + tempi[zonalindex]=tempiloc; //ion temperature + tempe[zonalindex]=tempeloc; //electron temperature + ui[zonalindex]=uiloc; //ion energy density + ue[zonalindex]=ueloc; //electron energy density + + #ifdef RELELECTRONS + gammabrk[zonalindex]=gbrkloc; + urelel[zonalindex]=urelelloc; + nrelel[zonalindex]=nrelelloc; + neth[zonalindex]=nethloc; + uratio_tot[zonalindex]=urelelloc/uint[zonalindex]; + uratio_th[zonalindex]=urelelloc/ueloc; + G0relel[zonalindex]=G0relelloc; + G0icrelel[zonalindex]=G0ic_relel_loc; + G0synrelel[zonalindex]=G0syn_relel_loc; + #endif + #endif + + ldouble temploc=calc_PEQ_Tfromurho(uint[zonalindex],rho[zonalindex],ix,iy,iz); + temp[zonalindex]=temploc; + entrlnT[zonalindex]=kB_over_mugas_mp*entr[zonalindex]/pp[RHO]; #ifdef CGSOUTPUT - rho[zonalindex]=rhoGU2CGS(rho[zonalindex]); - uint[zonalindex]=endenGU2CGS(uint[zonalindex]); - - #ifdef MAGNFIELD //ANDREW - bsq[zonalindex]=4.*M_PI*endenGU2CGS(bsq[zonalindex]); - #endif - - #ifdef PRINTVISCHEATINGTOSILO - vischeat[zonalindex]=endenGU2CGS(vischeat[zonalindex])*timeCGS2GU(1.); - vischeatnege[zonalindex]=endenGU2CGS(vischeatnege[zonalindex])*timeCGS2GU(1.); - vischeatnegi[zonalindex]=endenGU2CGS(vischeatnegi[zonalindex])*timeCGS2GU(1.); - #endif - #endif //CGSOUTPUT - - #ifdef PRINTCOULOMBTOSILO - coulomb[zonalindex]=calc_CoulombCoupling(pp,&geomout); - #endif - #ifdef PRINTGAMMATOSILO - gammag[zonalindex]=gamma; -#endif - - //default, but can be non-ortonormal VEL4 - lorentz[zonalindex]=fabs(vel[0])/sqrt(fabs(geomout.GG[0][0])); - u0[zonalindex]=fabs(vel[0]); - vx[zonalindex]=vel[1]; - vy[zonalindex]=vel[2]; - vz[zonalindex]=vel[3]; - - Edotx[zonalindex]=Tit[1]; - Edoty[zonalindex]=Tit[2]; - Edotz[zonalindex]=Tit[3]; - - //transform vel to cartesian - if (MYCOORDS==SCHWCOORDS || MYCOORDS==KSCOORDS || MYCOORDS==KERRCOORDS || MYCOORDS==SPHCOORDS || - MYCOORDS==MKS1COORDS || MYCOORDS==MKS2COORDS || MYCOORDS==MKS3COORDS || MYCOORDS==JETCOORDS || - MYCOORDS==MSPH1COORDS || MYCOORDS==MKER1COORDS) - { - vel[2]*=r; - vel[3]*=r*sin(th); - - vx[zonalindex] = sin(th)*cos(ph)*vel[1] - + cos(th)*cos(ph)*vel[2] - - sin(ph)*vel[3]; - - vy[zonalindex] = sin(th)*sin(ph)*vel[1] - + cos(th)*sin(ph)*vel[2] - + cos(ph)*vel[3]; - - vz[zonalindex] = cos(th)*vel[1] - - sin(th)*vel[2]; - - Tit[2]*=r; - Tit[3]*=r*sin(th); - - Edotx[zonalindex] = sin(th)*cos(ph)*Tit[1] - + cos(th)*cos(ph)*Tit[2] - - sin(ph)*Tit[3]; - - Edoty[zonalindex] = sin(th)*sin(ph)*Tit[1] - + cos(th)*sin(ph)*Tit[2] - + cos(ph)*Tit[3]; - - Edotz[zonalindex] = cos(th)*Tit[1] - - sin(th)*Tit[2]; - } - - - #ifdef MAGNFIELD - //magnetic field - Bx[zonalindex]=bcon[1]; - By[zonalindex]=bcon[2]; - Bz[zonalindex]=bcon[3]; - - #ifdef BATTERY - ldouble dBdtbat[4]; - estimate_Bgrowth_battery(ix,iy,iz,dBdtbat); - dBxdtbat[zonalindex]=dBdtbat[1]; - dBydtbat[zonalindex]=dBdtbat[2]; - dBzdtbat[zonalindex]=dBdtbat[3]; - #endif - - #ifdef MIMICDYNAMO - Bxdyn[zonalindex]=bcondyn[1]; - Bydyn[zonalindex]=bcondyn[2]; - Bzdyn[zonalindex]=bcondyn[3]; - #endif - - int iphimin,iphimax; - iphimin=0; - iphimax=ny-1; - #if defined(CORRECT_POLARAXIS) || defined(CORRECT_POLARAXIS_3D) - iphimin=NCCORRECTPOLAR; - #ifndef HALFTHETA - iphimax=ny-NCCORRECTPOLAR-1; - #endif - #endif - if(iy==iphimin) - { - phi[zonalindex]=geom.gdet*get_u(p,B1,ix,iy,iz)*get_size_x(iy,1)*2.*M_PI; - } - else if(iy>iphimin && iy<=iphimax) - { - imz=iz;imy=iy;imx=ix; -#ifdef PRINTXGC_RIGHT - imx=ix-NG; -#endif -#ifdef PRINTXGC_LEFT - imx=ix-NG; -#endif -#ifdef PRINTYGC_RIGHT - imy=iy-NG; -#endif - int idx=imz*(ny*nx) + (imy-1)*nx + imx; - phi[zonalindex]=phi[idx]+geom.gdet*get_u(p,B1,ix,iy,iz)*get_size_x(iy,1)*2.*M_PI; - } - - #ifdef MIMICDYNAMO - Avecdyn[zonalindex]=get_u(ptemp1,B3,ix,iy,iz); - if(iy==iphimin) - { - phidyn[zonalindex]=geom.gdet*get_u(pvecpot,1,ix,iy,iz)*get_size_x(iy,1)*2.*M_PI; - } - else if(iy>iphimin && iy<=iphimax) - { - imz=iz;imy=iy;imx=ix; -#ifdef PRINTXGC_RIGHT - imx=ix-NG; -#endif -#ifdef PRINTXGC_LEFT - imx=ix+NG; -#endif -#ifdef PRINTYGC_RIGHT - imy=iy-NG; -#endif - int idx=imz*(ny*nx) + (imy-1)*nx + imx; - phidyn[zonalindex]=phidyn[idx]+geom.gdet*get_u(pvecpot,1,ix,iy,iz)*get_size_x(iy,1)*2.*M_PI; - } - #endif - - #ifdef BATTERY - if(iy==iphimin) - { - phibat[zonalindex]=geom.gdet*dBdtbat[1]*get_size_x(iy,1)*2.*M_PI; - } - else if(iy>iphimin && iy<=iphimax) - { - imz=iz;imy=iy;imx=ix; -#ifdef PRINTXGC_RIGHT - imx=ix-NG; -#endif -#ifdef PRINTXGC_LEFT - imx=ix+NG; + rho[zonalindex]=rhoGU2CGS(rho[zonalindex]); + uint[zonalindex]=endenGU2CGS(uint[zonalindex]); + + #ifdef MAGNFIELD + bsq[zonalindex]=4.*M_PI*endenGU2CGS(bsq[zonalindex]); + #endif + + #ifdef EVOLVEELECTRONS + ui[zonalindex]=endenGU2CGS(ui[zonalindex]); + ue[zonalindex]=endenGU2CGS(ue[zonalindex]); + + #ifdef RELELECTRONS + urelel[zonalindex]=endenGU2CGS(urelel[zonalindex]); + nrelel[zonalindex]=numdensGU2CGS(nrelel[zonalindex]); + neth[zonalindex]=numdensGU2CGS(neth[zonalindex]); + G0relel[zonalindex]=endenGU2CGS(G0relel[zonalindex])*timeCGS2GU(1.); + G0icrelel[zonalindex]=endenGU2CGS(G0icrelel[zonalindex])*timeCGS2GU(1.); + G0synrelel[zonalindex]=endenGU2CGS(G0synrelel[zonalindex])*timeCGS2GU(1.); + #endif + #endif + + #ifdef PRINTVISCHEATINGTOSILO + vischeat[zonalindex] =endenGU2CGS(vischeat[zonalindex])*timeCGS2GU(1.); + vischeatnege[zonalindex]=endenGU2CGS(vischeatnege[zonalindex])*timeCGS2GU(1.); + vischeatnegi[zonalindex]=endenGU2CGS(vischeatnegi[zonalindex])*timeCGS2GU(1.); + #endif +#endif //CGSOUTPUT + +#ifdef PRINTCOULOMBTOSILO + coulomb[zonalindex]=calc_CoulombCoupling(pp,&geomout); #endif -#ifdef PRINTYGC_RIGHT - imy=iy-NG; + +#ifdef PRINTGAMMATOSILO + gammag[zonalindex]=gamma; +#endif + + // vector velocities and energy fluxes + vx[zonalindex]=vel[1]; + vy[zonalindex]=vel[2]; + vz[zonalindex]=vel[3]; + + Edotx[zonalindex]=Tit[1]; + Edoty[zonalindex]=Tit[2]; + Edotz[zonalindex]=Tit[3]; + + // transform velocity to cartesian + if (MYCOORDS==SCHWCOORDS || MYCOORDS==KSCOORDS || MYCOORDS==KERRCOORDS || MYCOORDS==SPHCOORDS || + MYCOORDS==MKS1COORDS || MYCOORDS==MKS2COORDS || MYCOORDS==MKS3COORDS || MYCOORDS==JETCOORDS || + MYCOORDS==MSPH1COORDS || MYCOORDS==MKER1COORDS) + { + vel[2]*=r; + vel[3]*=r*sin(th); + + vx[zonalindex] = sin(th)*cos(ph)*vel[1] + + cos(th)*cos(ph)*vel[2] + - sin(ph)*vel[3]; + + vy[zonalindex] = sin(th)*sin(ph)*vel[1] + + cos(th)*sin(ph)*vel[2] + + cos(ph)*vel[3]; + + vz[zonalindex] = cos(th)*vel[1] + - sin(th)*vel[2]; + + Tit[2]*=r; + Tit[3]*=r*sin(th); + + Edotx[zonalindex] = sin(th)*cos(ph)*Tit[1] + + cos(th)*cos(ph)*Tit[2] + - sin(ph)*Tit[3]; + + Edoty[zonalindex] = sin(th)*sin(ph)*Tit[1] + + cos(th)*sin(ph)*Tit[2] + + cos(ph)*Tit[3]; + + Edotz[zonalindex] = cos(th)*Tit[1] + - sin(th)*Tit[2]; + } + + // vector magnetic field quantities +#ifdef MAGNFIELD + //magnetic field + Bx[zonalindex]=bcon[1]; + By[zonalindex]=bcon[2]; + Bz[zonalindex]=bcon[3]; + + #ifdef BATTERY // battery field + ldouble dBdtbat[4]; + estimate_Bgrowth_battery(ix,iy,iz,dBdtbat); + dBxdtbat[zonalindex]=dBdtbat[1]; + dBydtbat[zonalindex]=dBdtbat[2]; + dBzdtbat[zonalindex]=dBdtbat[3]; + #endif + + #ifdef MIMICDYNAMO // dynamo field + ldouble ppdyn[NV]; int idyn; + PLOOP(idyn) ppdyn[idyn]=get_u(p,idyn,ix,iy,iz); + ppdyn[B1]=get_u(pvecpot,1,ix,iy,iz); + ppdyn[B2]=get_u(pvecpot,2,ix,iy,iz); + ppdyn[B3]=get_u(pvecpot,3,ix,iy,iz); + + #ifdef PRECOMPUTE_MY2OUT + trans_pmhd_coco_my2out(ppdyn, ppdyn, &geom, &geomout); + #else + trans_pmhd_coco(ppdyn, ppdyn, MYCOORDS,OUTCOORDS, xxvec,&geom,&geomout); + #endif + + ldouble bcondyn[4],bcovdyn[4]; + calc_bcon_prim(ppdyn,bcondyn,&geomout); + //indices_21(bcondyn,bcovdyn,geomout.gg); + + Bxdyn[zonalindex]=bcondyn[1]; + Bydyn[zonalindex]=bcondyn[2]; + Bzdyn[zonalindex]=bcondyn[3]; + #endif + + // calculate vector potential + //ANDREW changed phi defn to match harmpi + int iphimin,iphimax; + iphimin=0; + iphimax=ny-1; + + #if defined(CORRECT_POLARAXIS) || defined(CORRECT_POLARAXIS_3D) + //iphimin=NCCORRECTPOLAR; + #ifndef HALFTHETA + //iphimax=ny-NCCORRECTPOLAR-1; + #endif + #endif + if(iy==iphimin) + { + //phi[zonalindex]=geom.gdet*get_u(p,B1,ix,iy,iz)*get_size_x(iy,1)*2.*M_PI; + phi_accum[zonalindex]=geom.gdet*get_u(p,B1,ix,iy,iz)*get_size_x(iy,1);//*2.*M_PI; + phi[zonalindex] = phi_accum[zonalindex] - 0.5*geom.gdet*get_u(p,B1,ix,iy,iz)*get_size_x(iy,1); + } + else if(iy>iphimin && iy<=iphimax) + { + imz=iz;imy=iy;imx=ix; + #ifdef PRINTXGC_RIGHT + imx=ix-NG; + #endif + #ifdef PRINTXGC_LEFT + imx=ix-NG; + #endif + #ifdef PRINTYGC_RIGHT + imy=iy-NG; + #endif + int idx=imz*(ny*nx) + (imy-1)*nx + imx; + //phi[zonalindex]=phi[idx]+geom.gdet*get_u(p,B1,ix,iy,iz)*get_size_x(iy,1)*2.*M_PI; + phi_accum[zonalindex]=phi_accum[idx] + geom.gdet*get_u(p,B1,ix,iy,iz)*get_size_x(iy,1); //*2.*M_PI; + phi[zonalindex] = phi_accum[zonalindex] - 0.5*geom.gdet*get_u(p,B1,ix,iy,iz)*get_size_x(iy,1); + } + +#ifdef MIMICDYNAMO + Avecdyn[zonalindex]=get_u(ptemp1,B3,ix,iy,iz); + if(iy==iphimin) + { + phidyn[zonalindex]=geom.gdet*get_u(pvecpot,1,ix,iy,iz)*get_size_x(iy,1)*2.*M_PI; + } + else if(iy>iphimin && iy<=iphimax) + { + imz=iz;imy=iy;imx=ix; + #ifdef PRINTXGC_RIGHT + imx=ix-NG; + #endif + #ifdef PRINTXGC_LEFT + imx=ix+NG; + #endif + #ifdef PRINTYGC_RIGHT + imy=iy-NG; + #endif + int idx=imz*(ny*nx) + (imy-1)*nx + imx; + phidyn[zonalindex]=phidyn[idx]+geom.gdet*get_u(pvecpot,1,ix,iy,iz)*get_size_x(iy,1)*2.*M_PI; + } #endif - int idx=imz*(ny*nx) + (imy-1)*nx + imx; - phibat[zonalindex]=phibat[idx]+geom.gdet*dBdtbat[1]*get_size_x(iy,1)*2.*M_PI; - } - #endif - - //transform Bfield to cartesian - if (MYCOORDS==SCHWCOORDS || MYCOORDS==KSCOORDS || MYCOORDS==KERRCOORDS || MYCOORDS==SPHCOORDS || - MYCOORDS==MKS1COORDS || MYCOORDS==MKS2COORDS || MYCOORDS==MKS3COORDS || MYCOORDS==JETCOORDS || - MYCOORDS==MSPH1COORDS || MYCOORDS==MKER1COORDS) - { - bcon[2]*=r; - bcon[3]*=r*sin(th); - - Bx[zonalindex] = sin(th)*cos(ph)*bcon[1] - + cos(th)*cos(ph)*bcon[2] - - sin(ph)*bcon[3]; - - By[zonalindex] = sin(th)*sin(ph)*bcon[1] - + cos(th)*sin(ph)*bcon[2] - + cos(ph)*bcon[3]; - - Bz[zonalindex] = cos(th)*bcon[1] - - sin(th)*bcon[2]; #ifdef BATTERY - dBdtbat[2]*=r; - dBdtbat[3]*=r*sin(th); + if(iy==iphimin) + { + phibat[zonalindex]=geom.gdet*dBdtbat[1]*get_size_x(iy,1)*2.*M_PI; + } + else if(iy>iphimin && iy<=iphimax) + { + imz=iz;imy=iy;imx=ix; + #ifdef PRINTXGC_RIGHT + imx=ix-NG; + #endif + #ifdef PRINTXGC_LEFT + imx=ix+NG; + #endif + #ifdef PRINTYGC_RIGHT + imy=iy-NG; + #endif + int idx=imz*(ny*nx) + (imy-1)*nx + imx; + phibat[zonalindex]=phibat[idx]+geom.gdet*dBdtbat[1]*get_size_x(iy,1)*2.*M_PI; + } +#endif + + //transform Bfield to cartesian + if (MYCOORDS==SCHWCOORDS || MYCOORDS==KSCOORDS || MYCOORDS==KERRCOORDS || MYCOORDS==SPHCOORDS || + MYCOORDS==MKS1COORDS || MYCOORDS==MKS2COORDS || MYCOORDS==MKS3COORDS || MYCOORDS==JETCOORDS || + MYCOORDS==MSPH1COORDS || MYCOORDS==MKER1COORDS) + { + bcon[2]*=r; + bcon[3]*=r*sin(th); + + Bx[zonalindex] = sin(th)*cos(ph)*bcon[1] + + cos(th)*cos(ph)*bcon[2] + - sin(ph)*bcon[3]; + + By[zonalindex] = sin(th)*sin(ph)*bcon[1] + + cos(th)*sin(ph)*bcon[2] + + cos(ph)*bcon[3]; + + Bz[zonalindex] = cos(th)*bcon[1] + - sin(th)*bcon[2]; + +#ifdef BATTERY + dBdtbat[2]*=r; + dBdtbat[3]*=r*sin(th); - dBxdtbat[zonalindex] = sin(th)*cos(ph)*dBdtbat[1] - + cos(th)*cos(ph)*dBdtbat[2] - - sin(ph)*dBdtbat[3]; + dBxdtbat[zonalindex] = sin(th)*cos(ph)*dBdtbat[1] + + cos(th)*cos(ph)*dBdtbat[2] + - sin(ph)*dBdtbat[3]; - dBydtbat[zonalindex] = sin(th)*sin(ph)*dBdtbat[1] - + cos(th)*sin(ph)*dBdtbat[2] - + cos(ph)*dBdtbat[3]; + dBydtbat[zonalindex] = sin(th)*sin(ph)*dBdtbat[1] + + cos(th)*sin(ph)*dBdtbat[2] + + cos(ph)*dBdtbat[3]; - dBzdtbat[zonalindex] = cos(th)*dBdtbat[1] - - sin(th)*dBdtbat[2]; + dBzdtbat[zonalindex] = cos(th)*dBdtbat[1] + - sin(th)*dBdtbat[2]; #endif - #ifdef MIMICDYNAMO - bcondyn[2]*=r; - bcondyn[3]*=r*sin(th); +#ifdef MIMICDYNAMO + bcondyn[2]*=r; + bcondyn[3]*=r*sin(th); - Bxdyn[zonalindex] = sin(th)*cos(ph)*bcondyn[1] - + cos(th)*cos(ph)*bcondyn[2] - - sin(ph)*bcondyn[3]; + Bxdyn[zonalindex] = sin(th)*cos(ph)*bcondyn[1] + + cos(th)*cos(ph)*bcondyn[2] + - sin(ph)*bcondyn[3]; - Bydyn[zonalindex] = sin(th)*sin(ph)*bcondyn[1] - + cos(th)*sin(ph)*bcondyn[2] - + cos(ph)*bcondyn[3]; + Bydyn[zonalindex] = sin(th)*sin(ph)*bcondyn[1] + + cos(th)*sin(ph)*bcondyn[2] + + cos(ph)*bcondyn[3]; - Bzdyn[zonalindex] = cos(th)*bcondyn[1] - - sin(th)*bcondyn[2]; - #endif - } + Bzdyn[zonalindex] = cos(th)*bcondyn[1] + - sin(th)*bcondyn[2]; +#endif + } #endif //MAGNFIELD #ifdef RADIATION - ldouble Rtt,ehat,ugas[4],urad[4],rvel[4],Rij[4][4],Rij22[4][4],Giff[4],tradloc,tradlteloc; - struct opacities opac; - ldouble tauabsloc = vcon[0]*calc_kappa(pp,&geomout,&opac); - ldouble tauscaloc = vcon[0]*calc_kappaes(pp,&geomout); - ldouble taueffloc = sqrt(tauabsloc*(tauabsloc+tauscaloc)); - - if(doingavg==0) //from snapshot - { - calc_ff_Rtt(pp,&Rtt,ugas,&geomout); - ehat=-Rtt; - calc_Rij(pp,&geomout,Rij); //calculates R^munu in OUTCOORDS - - //four fource - - calc_Gi(pp,&geomout,Giff,0.0,0,0); //ANDREW fluid frame - - #ifdef RADFLUXFFINOUTPUT - boost22_lab2ff(Rij,Rij,pp,geomout.gg,geomout.GG); - #endif - - indices_2221(Rij,Rij,geomout.gg); - //calculating radiation temperature - tradlteloc=calc_LTE_TfromE(ehat); - tradloc=calc_LTE_TfromE(ehat); - -#ifdef EVOLVEPHOTONNUMBER - tradloc=calc_ncompt_Thatrad(pp,&geomout,ehat); -#endif - } - else //using avg - { - ehat=get_uavg(pavg,AVGEHAT,ix,iy,iz); - tradloc=calc_LTE_TfromE(ehat); - tradlteloc=calc_LTE_TfromE(ehat); - for(j=0;j<4;j++) - Giff[j]=get_uavg(pavg,AVGGHAT(j),ix,iy,iz); - - //ANDREW if Gi0 accidentally saved in lab frame: - #ifdef SIMOUTPUT_GILAB2FF - boost2_lab2ff(Giff,Giff,pp,geomout.gg,geomout.GG); //ANDREW avg Gff already in OUTCOORDS - #endif -#ifdef EVOLVEPHOTONNUMBER - tradloc=calc_ncompt_Thatrad_fromEN(ehat,get_uavg(pavg,AVGNFHAT,ix,iy,iz)); -#endif - for(i=0;i<4;i++) - for(j=0;j<4;j++) - Rij[i][j]=get_uavg(pavg,AVGRIJ(i,j),ix,iy,iz); - } - - //Bernoulli number - muBe[zonalindex]+=-Rij[1][0]/rhouconr; - Be[zonalindex]+=-Rij[0][0]/rhoucont; - - indices_2122(Rij,Rij22,geomout.GG); - - //correcting rad-velocities basing on - int radcorr; - - p2u(pp,uu,&geomout); - - uu[EE0]=gdetu*Rij[0][0]; - uu[FX0]=gdetu*Rij[0][1]; - uu[FY0]=gdetu*Rij[0][2]; - uu[FZ0]=gdetu*Rij[0][3]; - - //print_conserved(uu); - u2p_rad(uu,pp,&geomout,&radcorr); - - #ifdef COMPTONIZATION - //Calculate thermal comptonization - ldouble Gic_tmp[4]; - ldouble ucon_ff[4]; - ucon_ff[1]=ucon_ff[2]=ucon_ff[3]=0.; - ucon_ff[0]=1.; - ldouble kappaes=calc_kappaes(pp, &geomout); - calc_Compt_Gi(pp, &geomout, Gic_tmp, ehat, tempeloc, kappaes, ucon_ff); - G0ic_th_loc = fabs(Gic_tmp[0]); - - G0icth[zonalindex]=G0ic_th_loc; - #ifdef CGSOUTPUT - G0icth[zonalindex]=endenGU2CGS(G0icth[zonalindex]) * timeCGS2GU(1.); - #endif - #endif - - Ehat[zonalindex]=ehat; - Gtff[zonalindex]=-Giff[0]; - Erad[zonalindex]=Rij22[0][0]; - - #ifdef EVOLVEPHOTONNUMBER - Nph[zonalindex]=pp[NF0]; - #endif - - #ifdef CGSOUTPUT - Ehat[zonalindex]=endenGU2CGS(Ehat[zonalindex]); - Gtff[zonalindex]=endenGU2CGS(Gtff[zonalindex])*timeCGS2GU(1.); - Erad[zonalindex]=endenGU2CGS(Erad[zonalindex]); - #ifdef EVOLVEPHOTONNUMBER - Nph[zonalindex]=numdensGU2CGS(Nph[zonalindex]); - #endif - #endif - - trad[zonalindex]=tradloc; - tradlte[zonalindex]=tradlteloc; - entrrad[zonalindex]=(4./3.)*A_RAD*tradloc*tradloc*tradloc/pp[RHO]; - - Fx[zonalindex]=Rij22[1][0]; - Fy[zonalindex]=Rij22[2][0]; - Fz[zonalindex]=Rij22[3][0]; - - urad[1]=pp[FX0]; - urad[2]=pp[FY0]; - urad[3]=pp[FZ0]; - conv_vels(urad,urad,VELPRIM,VEL4,geomout.gg,geomout.GG); - - uradx[zonalindex]=urad[1]; - urady[zonalindex]=urad[2]; - uradz[zonalindex]=urad[3]; - - //Angular integration (tau_theta) - if(iy==0) - { - tausca[zonalindex]=tauscaloc*dxph[1]; - tauabs[zonalindex]=tauabsloc*dxph[1]; - taueff[zonalindex]=taueffloc*dxph[1]; - } - else - { - imz=iz; - imy=iy; - imx=ix; -#ifdef PRINTXGC_RIGHT - imx=ix-NG; -#endif -#ifdef PRINTXGC_LEFT - imx=ix+NG; -#endif -#ifdef PRINTYGC_RIGHT - imy=iy-NG; -#endif - int idx=imz*(ny*nx) + (imy-1)*nx + imx; - if(iy<=NY/2) //proper integration only in the upper half - { - tausca[zonalindex]=tausca[idx]+tauscaloc*dxph[1]; - tauabs[zonalindex]=tauabs[idx]+tauabsloc*dxph[1]; - taueff[zonalindex]=taueff[idx]+taueffloc*dxph[1]; - } - else - { - idx=imz*(ny*nx) + (NY/2-1)*nx + imx; - tausca[zonalindex]=tausca[idx]; - tauabs[zonalindex]=tauabs[idx]; - taueff[zonalindex]=taueff[idx]; - } - } - - - //transform rad flux to cartesian - if (MYCOORDS==SCHWCOORDS || MYCOORDS==KSCOORDS || MYCOORDS==KERRCOORDS || MYCOORDS==SPHCOORDS || - MYCOORDS==MKS1COORDS || MYCOORDS==MKS2COORDS || MYCOORDS==MKS3COORDS || MYCOORDS==JETCOORDS || - MYCOORDS==MSPH1COORDS || MYCOORDS==MKER1COORDS) - { - Rij22[2][0]*=sqrt(geomout.gg[2][2]); - Rij22[3][0]*=sqrt(geomout.gg[3][3]); - - //Rij22? - - Fx[zonalindex] = sin(th)*cos(ph)*Rij22[1][0] - + cos(th)*cos(ph)*Rij22[2][0] - - sin(ph)*Rij22[3][0]; - - Fy[zonalindex] = sin(th)*sin(ph)*Rij22[1][0] - + cos(th)*sin(ph)*Rij22[2][0] - + cos(ph)*Rij22[3][0]; - - Fz[zonalindex] = cos(th)*Rij22[1][0] - - sin(th)*Rij22[2][0]; - - urad[2]*=r; - urad[3]*=r*sin(th); - - uradx[zonalindex] = sin(th)*cos(ph)*urad[1] - + cos(th)*cos(ph)*urad[2] - - sin(ph)*urad[3]; - - urady[zonalindex] = sin(th)*sin(ph)*urad[1] - + cos(th)*sin(ph)*urad[2] - + cos(ph)*urad[3]; - - uradz[zonalindex] = cos(th)*urad[1] - - sin(th)*urad[2]; - } - -#endif //RADIATION - - - + ldouble Rtt,ehat,ugas[4],urad[4],rvel[4],Rij[4][4],Rij22[4][4],Giff[4],tradloc,tradlteloc; + struct opacities opac; + ldouble tauabsloc = vcon[0]*calc_kappa(pp,&geomout,&opac); + ldouble tauscaloc = vcon[0]*calc_kappaes(pp,&geomout); + ldouble taueffloc = sqrt(tauabsloc*(tauabsloc+tauscaloc)); + + if(doingavg==0) //from snapshot + { + calc_ff_Rtt(pp,&Rtt,ugas,&geomout); + ehat=-Rtt; + calc_Rij(pp,&geomout,Rij); //calculates R^munu in OUTCOORDS + + //four fource + calc_Gi(pp,&geomout,Giff,0.0,0,0); // fluid frame + + #ifdef RADFLUXFFINOUTPUT + boost22_lab2ff(Rij,Rij,pp,geomout.gg,geomout.GG); + #endif + indices_2221(Rij,Rij,geomout.gg); + + //calculating radiation temperature + tradlteloc=calc_LTE_TfromE(ehat); + tradloc=calc_LTE_TfromE(ehat); + + #ifdef EVOLVEPHOTONNUMBER + tradloc=calc_ncompt_Thatrad(pp,&geomout,ehat); + #endif + } + else //using avg + { + ehat=get_uavg(pavg,AVGEHAT,ix,iy,iz); + tradloc=calc_LTE_TfromE(ehat); + tradlteloc=calc_LTE_TfromE(ehat); + for(j=0;j<4;j++) + Giff[j]=get_uavg(pavg,AVGGHAT(j),ix,iy,iz); + + //ANDREW if Gi0 accidentally saved in lab frame: + #ifdef SIMOUTPUT_GILAB2FF + boost2_lab2ff(Giff,Giff,pp,geomout.gg,geomout.GG); //ANDREW avg Gff already in OUTCOORDS + #endif + #ifdef EVOLVEPHOTONNUMBER + tradloc=calc_ncompt_Thatrad_fromEN(ehat,get_uavg(pavg,AVGNFHAT,ix,iy,iz)); + #endif + for(i=0;i<4;i++) + for(j=0;j<4;j++) + Rij[i][j]=get_uavg(pavg,AVGRIJ(i,j),ix,iy,iz); + } + indices_2122(Rij,Rij22,geomout.GG); + + // Add radiation to Bernoulli number + muBe[zonalindex]+=-Rij[1][0]/rhouconr; + Be[zonalindex]+=-Rij[0][0]/rhoucont; + + //correcting rad-velocities basing on + int radcorr; + + p2u(pp,uu,&geomout); + + uu[EE0]=gdetu*Rij[0][0]; + uu[FX0]=gdetu*Rij[0][1]; + uu[FY0]=gdetu*Rij[0][2]; + uu[FZ0]=gdetu*Rij[0][3]; + + //print_conserved(uu); + u2p_rad(uu,pp,&geomout,&radcorr); + + #ifdef COMPTONIZATION + //Calculate thermal comptonization + ldouble Gic_tmp[4]; + ldouble ucon_ff[4]; + ucon_ff[1]=ucon_ff[2]=ucon_ff[3]=0.; + ucon_ff[0]=1.; + ldouble kappaes=calc_kappaes(pp, &geomout); + calc_Compt_Gi(pp, &geomout, Gic_tmp, ehat, tempeloc, kappaes, ucon_ff); + G0ic_th_loc = fabs(Gic_tmp[0]); + + G0icth[zonalindex]=G0ic_th_loc; + #endif + + Ehat[zonalindex]=ehat; + Gtff[zonalindex]=-Giff[0]; + Erad[zonalindex]=Rij22[0][0]; + + #ifdef EVOLVEPHOTONNUMBER + Nph[zonalindex]=pp[NF0]; + #endif + + #ifdef CGSOUTPUT + Ehat[zonalindex]=endenGU2CGS(Ehat[zonalindex]); + Gtff[zonalindex]=endenGU2CGS(Gtff[zonalindex])*timeCGS2GU(1.); + Erad[zonalindex]=endenGU2CGS(Erad[zonalindex]); + #ifdef EVOLVEPHOTONNUMBER + Nph[zonalindex]=numdensGU2CGS(Nph[zonalindex]); + #endif + #ifdef COMPTONIZATION + G0icth[zonalindex]=endenGU2CGS(G0icth[zonalindex]) * timeCGS2GU(1.); + #endif + #endif + + trad[zonalindex]=tradloc; + tradlte[zonalindex]=tradlteloc; + entrrad[zonalindex]=(4./3.)*A_RAD*tradloc*tradloc*tradloc/pp[RHO]; + + Fx[zonalindex]=Rij22[1][0]; + Fy[zonalindex]=Rij22[2][0]; + Fz[zonalindex]=Rij22[3][0]; + + urad[1]=pp[FX0]; + urad[2]=pp[FY0]; + urad[3]=pp[FZ0]; + conv_vels(urad,urad,VELPRIM,VEL4,geomout.gg,geomout.GG); + + uradx[zonalindex]=urad[1]; + urady[zonalindex]=urad[2]; + uradz[zonalindex]=urad[3]; + + //Angular integration (tau_theta) + if(iy==0) + { + tausca[zonalindex]=tauscaloc*dxph[1]; + tauabs[zonalindex]=tauabsloc*dxph[1]; + taueff[zonalindex]=taueffloc*dxph[1]; + } + else + { + imz=iz; + imy=iy; + imx=ix; + #ifdef PRINTXGC_RIGHT + imx=ix-NG; + #endif + #ifdef PRINTXGC_LEFT + imx=ix+NG; + #endif + #ifdef PRINTYGC_RIGHT + imy=iy-NG; + #endif + int idx=imz*(ny*nx) + (imy-1)*nx + imx; + if(iy<=NY/2) //proper integration only in the upper half + { + tausca[zonalindex]=tausca[idx]+tauscaloc*dxph[1]; + tauabs[zonalindex]=tauabs[idx]+tauabsloc*dxph[1]; + taueff[zonalindex]=taueff[idx]+taueffloc*dxph[1]; + } + else + { + idx=imz*(ny*nx) + (NY/2-1)*nx + imx; + tausca[zonalindex]=tausca[idx]; + tauabs[zonalindex]=tauabs[idx]; + taueff[zonalindex]=taueff[idx]; } } + + //transform rad flux to cartesian + if (MYCOORDS==SCHWCOORDS || MYCOORDS==KSCOORDS || MYCOORDS==KERRCOORDS || MYCOORDS==SPHCOORDS || + MYCOORDS==MKS1COORDS || MYCOORDS==MKS2COORDS || MYCOORDS==MKS3COORDS || MYCOORDS==JETCOORDS || + MYCOORDS==MSPH1COORDS || MYCOORDS==MKER1COORDS) + { + Rij22[2][0]*=sqrt(geomout.gg[2][2]); + Rij22[3][0]*=sqrt(geomout.gg[3][3]); + + Fx[zonalindex] = sin(th)*cos(ph)*Rij22[1][0] + + cos(th)*cos(ph)*Rij22[2][0] + - sin(ph)*Rij22[3][0]; + + Fy[zonalindex] = sin(th)*sin(ph)*Rij22[1][0] + + cos(th)*sin(ph)*Rij22[2][0] + + cos(ph)*Rij22[3][0]; + + Fz[zonalindex] = cos(th)*Rij22[1][0] + - sin(th)*Rij22[2][0]; + + urad[2]*=r; + urad[3]*=r*sin(th); + + uradx[zonalindex] = sin(th)*cos(ph)*urad[1] + + cos(th)*cos(ph)*urad[2] + - sin(ph)*urad[3]; + + urady[zonalindex] = sin(th)*sin(ph)*urad[1] + + cos(th)*sin(ph)*urad[2] + + cos(ph)*urad[3]; + + uradz[zonalindex] = cos(th)*urad[1] + - sin(th)*urad[2]; + } + +#endif //RADIATION + } } + } //Loop for radially integrated optical depth @@ -1263,223 +1221,204 @@ int fprint_silofile(ldouble time, int num, char* folder, char* prefix) #else for(iz=0;iz NG-1; ix--) + { +#if defined(PRINTXGC_RIGHT) + for(ix = nx+NG-1; ix > NG-1; ix--) #elif defined(PRINTXGC_LEFT) - for(ix = nx-NG-1;ix > -NG-1; ix--) + for(ix = nx-NG-1;ix > -NG-1; ix--) #else - for(ix=(nx-1);ix>-1;ix--) + for(ix=(nx-1);ix>-1;ix--) #endif - { + { #ifdef PRINTYGC_RIGHT for(iy=NG;iyix,geom->iy,geom->iz,&gix,&giy,&giz); - if(verbose) printf("%4d > %4d %4d %4d > NEGUU > neg uu[0] - requesting fixup\n",PROCID,gix,giy,giz); - pp[0]=RHOFLOOR; //used when not fixing up - uu[0]=RHOFLOOR*gdetu; - ret=-2; //to request fixup - //ANDREW -- but ret=-1 if energy inversion failes but entropy inversion does not! - //ANDREW -- do we always want a fixup if we have negative uu[0] ? - u2pret=-1; // indicates that inversion is needed + 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 + uu[RHO]=RHOFLOOR*gdetu; + ret=-2; // to request fixup in all cases #ifndef SWAPPAPC - global_int_slot[GLOBALINTSLOT_NTOTALMHDFIXUPS]++; //but count as fixup + global_int_slot[GLOBALINTSLOT_NTOTALMHDFIXUPS]++; // count as fixup #endif } - if(u2pret!=0) // u2pret=-1 at this stage, so this is always satisfied - { #ifdef ENFORCEENTROPY - u2pret=-1; //skip hot energy-conserving inversion and go to entropy inversion + u2pret = -1; //skip hot energy-conserving inversion and go to entropy inversion #else - u2pret = u2p_solver(uu,pp,ggg,U2P_HOT,0); // invert using the hot energy equation + u2pret = u2p_solver(uu,pp,ggg,U2P_HOT,0); // invert using the hot energy equation #endif //ENFORCEENTROPY - } - + if(ALLOWENTROPYU2P) // Inversion with entropy equation -- on by default (see choices.h) { if(u2pret<0) // true if energy equation failed, or if energy equation was not required (because ENFORCEENTROPY is defined) @@ -220,15 +205,12 @@ 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[1]); + printf("u2p_entr >>> %d %d <<< %d >>> %e > %e\n",geom->ix + TOI, geom->iy + TOJ,u2pret,u0,pp[UU]); } //************************************ - //entropy solver - conserving entropy - if(ppentr[RHO]<0.) //if not yet calculated - { - u2pret=u2p_solver(uu,pp,ggg,U2P_ENTROPY,0); // invert using entropy equation - } + //entropy solver - invert using entropy equation + u2pret = u2p_solver(uu,pp,ggg,U2P_ENTROPY,0); if(u2pret<0) { @@ -236,18 +218,17 @@ 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[0],uu[1],uu[5],pp[0],pp[1],geom->ix,geom->iy,geom->iz); + 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); } } // if(u2pret<0) // second time -- entropy eqn } // if(u2pret<0) // first time -- energy eqn - } // if(ALLOWENTROPYU2P) + } // if(ALLOWENTROPYU2P) if(u2pret<0) // entropy equation also failed { - - //leaving primitives unchanged - should not happen - if(verbose>1 || 1) + //leaving primitives unchanged + if(verbose > 1 || 1) { printf("%4d > %4d %4d %4d > MHDU2PFAIL > u2p prim. unchanged > %d \n",PROCID,geom->ix+TOI,geom->iy+TOJ,geom->iz+TOK,u2pret); } @@ -256,9 +237,9 @@ u2p(ldouble *uu0, ldouble *pp, void *ggg, int corrected[3], int fixups[2], int t pp[u2pret]=ppbak[u2pret]; } - if(ret<0) //to update conserved + if(ret<0) // update conserved hdcorr=1; - if(ret<-1) //request fixup when entropy failed + if(ret<-1) // request fixup when entropy failed fixups[0]=1; else fixups[0]=0; @@ -345,24 +326,30 @@ check_floors_mhd(ldouble *pp, int whichvel,void *ggg) ldouble uu[NV],pporg[NV]; p2u_mhd(pp,uu,ggg); - + + ldouble pgamma=GAMMA; + #ifdef CONSISTENTGAMMA + pgamma=pick_gammagas(geom->ix,geom->iy,geom->iz); + #endif + //********************************************************************** //rho too small - if(pp[0]ix+TOI,geom->iy+TOJ,geom->iz+TOK,geom->ix,geom->iy,geom->iz,pp[0],TI,TJ,TK); - pp[0]=RHOFLOOR; + if(verbose ) printf("hd_floors CASE 1 at %d %d %d | %d %d %d (%e) | tijk: %d %d %d\n",geom->ix+TOI,geom->iy+TOJ,geom->iz+TOK,geom->ix,geom->iy,geom->iz,pp[RHO],TI,TJ,TK); + pp[RHO]=RHOFLOOR; ret=-1; } //********************************************************************** - //rho too small, use floor scaling as r^-3/2 + //rho too small, + // floor scaling as r^-3/2 #ifdef RHOFLOOR_BH ldouble xxBL[4]; #ifdef PRECOMPUTE_MY2OUT @@ -373,22 +360,23 @@ check_floors_mhd(ldouble *pp, int whichvel,void *ggg) ldouble rr = xxBL[1] / rhorizonBL; ldouble rhofloor = RHOFLOOR_BH_NORM / sqrt(rr*rr*rr); - if(pp[0]ix+TOI,geom->iy+TOJ,pp[0]); - pp[0]=rhofloor; + if(verbose ) printf("hd_floors BH CASE 1 at %d %d (%e)\n",geom->ix+TOI,geom->iy+TOJ,pp[RHO]); + 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.; @@ -401,16 +389,16 @@ check_floors_mhd(ldouble *pp, int whichvel,void *ggg) ldouble rr = xxBL[1] / rout; ldouble rhofloor = RHOATMMIN / sqrt(rr*rr*rr); ldouble uintfloor = UINTATMMIN / sqrt(rr*rr*rr*rr*rr); - if(pp[0] < rhofloor || pp[1] < uintfloor) + if(pp[RHO] < rhofloor || pp[UU] < uintfloor) { for(iv = 0; iv < NVMHD; iv++) { pporg[iv] = pp[iv]; } - if(verbose ) printf("hd_floors BH CASE INIT at %d %d (%e)\n",geom->ix+TOI, geom->iy+TOJ, pp[0]); - pp[0] = rhofloor; - pp[1] = uintfloor; + if(verbose ) printf("hd_floors BH CASE INIT at %d %d (%e)\n",geom->ix+TOI, geom->iy+TOJ, pp[RHO]); + pp[RHO] = rhofloor; + pp[UU] = uintfloor; ret=-1; } @@ -418,15 +406,15 @@ check_floors_mhd(ldouble *pp, int whichvel,void *ggg) //********************************************************************** //too cold - if(pp[1]ix+TOI,geom->iy+TOJ,geom->iz+TOK,geom->ix,geom->iy,geom->iz,pp[0],pp[1],TI,TJ,TK);}//getchar();} - pp[1]=UURHORATIOMIN*pp[0]; //increasing uint + 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[UU]; //increasing uint ret=-1; @@ -434,124 +422,232 @@ check_floors_mhd(ldouble *pp, int whichvel,void *ggg) //********************************************************************** //too hot - if(pp[1]>UURHORATIOMAX*pp[0]) + if(pp[UU]>UURHORATIOMAX*pp[RHO]) { for(iv=0;ivix+TOI,geom->iy+TOJ,geom->iz,pp[0],pp[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]); } //********************************************************************** //too magnetized #ifdef MAGNFIELD - ldouble ucond[4],ucovd[4]; - ldouble bcond[4],bcovd[4],bsq,magpre; + ldouble ucon[4],ucov[4]; + ldouble bcon[4],bcov[4],bsq,magpre; ldouble etacon[4],etarel[4]; + ldouble f=1.; + ldouble fuu=1.; + int rhofloored=0; + int uufloored=0; + for(iv=1;iv<4;iv++) - ucond[iv]=pp[1+iv]; - calc_ucon_ucov_from_prims(pp, geom, ucond, ucovd); - calc_bcon_bcov_bsq_from_4vel(pp, ucond, ucovd, geom, bcond, bcovd, &bsq); - //magpre = 0.5 * bsq; - magpre = bsq; // ANDREW: B2RHORATIOMAX and B2UURATIOMAX make more sense with bsq, not .5bsq - - calc_normalobs_ncon(GG, geom->alpha, etacon); - conv_vels_ut(etacon,etarel,VEL4,VELPRIM,gg,GG); + 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; + //magpre = bsq; // ANDREW: B2RHORATIOMAX and B2UURATIOMAX make more sense with bsq, not 0.5 bsq + + // 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); + f=bsq/(B2RHORATIOMAX*pp[RHO]); //increase factor, f>1 + rhofloored=1; + } - if(magpre>B2RHORATIOMAX*pp[RHO]) + // 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[RHO],magpre); - ldouble f=magpre/(B2RHORATIOMAX*pp[RHO]); + 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 for(iv=0;ivbetasqmax) + { + betasq=betasqmax; + printf("floored parallel velocity\n"); + } + gammapar = 1./sqrt(1.-betasq); + + // drift frame velocity + for(iv=0;iv<4;iv++) + ucondr[iv] = gammapar*(ucon[iv] + betapar*bcon[iv]); + + // magnetic field + Bcon[0]=0.; + Bcon[1]=pp[B1]; + Bcon[2]=pp[B2]; + Bcon[3]=pp[B3]; + indices_21(Bcon,Bcov,gg); + Bsq = dot(Bcon,Bcov); + Bmag = sqrt(Bsq); + + // conserved parallel momentum + udotB = dot(ucon,Bcov); + QdotB = udotB*wold*ucon[0]; + + // get new parallel velocity + xx = 2.*QdotB/(Bmag*wnew*ucondr[0]); + vpar = xx / (ucondr[0]*(1.+sqrt(1.+xx*xx))); + + // new three velocity + vcon[0]=1.; + for(iv=1;iv<4;iv++) + { + vcon[iv]=vpar*Bcon[iv]/(Bmag) + ucondr[iv]/ucondr[0]; + } + + // to VELPRIM + conv_vels(vcon,ucont,VEL3,VELPRIM,gg,GG); + //conv_vels(uconout,uconout,VEL4,VELPRIM,gg,GG); + + pp[VX] = ucont[1]; + pp[VY] = ucont[2]; + pp[VZ] = ucont[3]; + + // check + calc_ucon_ucov_from_prims(pp, geom, ucon, ucov); + calc_bcon_bcov_bsq_from_4vel(pp, ucon, ucov, geom, bcon, bcov, &bsq); - 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 + + ldouble dpp[NV],duu[NV]; + ldouble drho=pp[RHO]*(f-1.); + ldouble dugas=pp[UU]*(fuu-1.); + + for(iv=0;ivalpha, etacon); + conv_vels_ut(etacon,etarel,VEL4,VELPRIM,gg,GG); + + dpp[RHO]=drho; + //ANDREW this is a change, used to keep UU fixed here + dpp[UU]=dugas; + //dpp[UU]=0.; + dpp[VX] = etarel[1]; + dpp[VY] = etarel[2]; + dpp[VZ] = etarel[3]; + dpp[ENTR] = 0.; + dpp[B1] = dpp[B2] = dpp[B3] = 0.; + + // delta conserved in ZAMO + p2u_mhd(dpp,duu,geom); - for(iv=0;ivix+TOI>5) //report only outside horizon -#endif - printf("u2p failed after imposing bsq over rho floors at %d %d %d with gamma=%f\n",geom->ix+TOI,geom->iy+TOJ,geom->iz+TOK,get_u_scalar(gammagas,geom->ix,geom->iy,geom->iz)); + 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); #ifdef B2RHOFLOOR_BACKUP_FFFRAME - // Backup bsq/rho floor -- if zamo frame fails, do fluid frame instead of crashing - for(iv=0;ivix,geom->iy,geom->iz); ldouble ne=calc_thermal_ne(pporg); //thermal only - ldouble pe=K_BOLTZ*ne*Te; - ldouble gammae=GAMMAE; -#ifdef CONSISTENTGAMMA -#ifndef FIXEDGAMMASPECIES - gammae=calc_gammaintfromtemp(Te,ELECTRONS); -#endif -#endif - ldouble ue=pe/(gammae-1.); ldouble ni=pporg[RHO]/MU_I/M_PROTON; + + ldouble pe=K_BOLTZ*ne*Te; ldouble pi=K_BOLTZ*ni*Ti; + + ldouble gammae=GAMMAE; ldouble gammai=GAMMAI; -#ifdef CONSISTENTGAMMA -#ifndef FIXEDGAMMASPECIES + #ifdef CONSISTENTGAMMA + #ifndef FIXEDGAMMASPECIES + gammae=calc_gammaintfromtemp(Te,ELECTRONS); gammai=calc_gammaintfromtemp(Ti,IONS); -#endif -#endif + #endif + #endif + + ldouble ue=pe/(gammae-1.); ldouble ui=pi/(gammai-1.); //calculate new entropy using pp[] @@ -562,55 +658,45 @@ check_floors_mhd(ldouble *pp, int whichvel,void *ggg) mass = M_ELECTR; gamma = GAMMAE; n = calc_thermal_ne(pp); -#ifdef CONSISTENTGAMMA + #ifdef CONSISTENTGAMMA Tenew=solve_Teifromnmu(n, mass, ue,ELECTRONS); //solves in parallel for gamma and temperature theta=K_BOLTZ*Tenew/mass; -#ifndef FIXEDGAMMASPECIES + #ifndef FIXEDGAMMASPECIES gamma=calc_gammaintfromtheta(theta); //the same gamma as just solved -#endif -#endif + #endif + #endif pnew=(ue)*(gamma-1.); Tenew=pnew/K_BOLTZ/n; - - ldouble rhoe=n*MU_E*M_PROTON; - Senew=calc_SefromrhoT(rhoe,Tenew,ELECTRONS); - + Senew=calc_SefromrhoT(n*MU_E*M_PROTON,Tenew,ELECTRONS); pp[ENTRE]=Senew; //ions mass = M_PROTON; gamma = GAMMAI; - n = calc_thermal_ne(pp); -#ifdef CONSISTENTGAMMA + 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 + #ifndef FIXEDGAMMASPECIES gamma=calc_gammaintfromtheta(theta); //the same gamma as just solved -#endif -#endif - + #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(magpre>B2RHORATIOMAX*pp[RHO]) + } //if(bsq>B2RHORATIOMAX*pp[RHO]) or (bsq>B2UURATIOMAX*pp[UU]) - //independent check on ugas vs - if(magpre>B2UURATIOMAX*pp[UU]) - { - //if(verbose) printf("mag_floors CASE 3 at (%d,%d,%d): %e %e\n",geom->ix+TOI,geom->iy+TOJ,geom->iz,pp[UU],magpre); - pp[UU]*=magpre/(B2UURATIOMAX*pp[UU]); - ret=-1; - } #endif //MAGNFIELD //********************************************************************** //too fast + //TODO: implement checks for other VELPRIM if(VELPRIM==VELR) { ldouble qsq=0.; @@ -639,7 +725,7 @@ check_floors_mhd(ldouble *pp, int whichvel,void *ggg) } } - //TODO: implement checks for other VELPRIM + //********************************************************************** //Species temperature floors/ceilings @@ -758,7 +844,7 @@ check_floors_mhd(ldouble *pp, int whichvel,void *ggg) //ANDREW TODO do we want this? Is this inconsistent with keeping entropy as a backup until the very end of time step? //updates entropy after floor corrections if(ret<0) - pp[5]=calc_Sfromu(pp[RHO],pp[UU],geom->ix,geom->iy,geom->iz); + pp[ENTR]=calc_Sfromu(pp[RHO],pp[UU],geom->ix,geom->iy,geom->iz); return ret; } @@ -1001,10 +1087,15 @@ f_u2p_entropy(ldouble Wp, ldouble* cons, ldouble *f, ldouble *df, ldouble *err,l int u2p_solver(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose) { - #ifdef NONRELMHD + +#ifdef NONRELMHD return u2p_solver_nonrel(uu,pp,ggg,Etype,verbose); - #endif +#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; @@ -1105,7 +1196,6 @@ u2p_solver_nonrel(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose) //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; @@ -1118,7 +1208,6 @@ u2p_solver_nonrel(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose) int ib; for(ib=0;ibalpha; //D - D = uu[0] * gdetu_inv * alpha; //uu[0]=gdetu rho ut + D = uu[RHO] * gdetu_inv * alpha; // gdetu rho ut //conserved entropy "S u^t" - Sc = uu[5] * gdetu_inv * alpha; + Sc = uu[ENTR] * gdetu_inv * alpha; //Q_mu=alpha T^t_mu - Qcov[0] = (uu[1] * gdetu_inv - uu[0] * gdetu_inv) * alpha; - Qcov[1] = uu[2] * gdetu_inv * alpha; - Qcov[2] = uu[3] * gdetu_inv * alpha; - Qcov[3] = uu[4] * gdetu_inv * alpha; + 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[1] * gdetu_inv *alpha; - Qcovp[1] = uu[2] * gdetu_inv *alpha; - Qcovp[2] = uu[3] * gdetu_inv *alpha; - Qcovp[3] = uu[4] * gdetu_inv *alpha; + 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); @@ -1693,9 +1782,6 @@ u2p_solver_W(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose) Qn = Qcon[0] * ncov[0]; //j^mu_nu = delta^mu_nu +n^mu n_nu -#ifdef APPLY_OMP_SIMD - //#pragma omp simd -#endif for(i=0;i<4;i++) for(j=0;j<4;j++) jmunu[i][j] = delta(i,j) + ncon[i]*ncov[j]; @@ -1704,9 +1790,6 @@ u2p_solver_W(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose) for(i=0;i<4;i++) { Qtcon[i]=0.; -#ifdef APPLY_OMP_SIMD - //#pragma omp simd -#endif for(j=0;j<4;j++) Qtcon[i]+=jmunu[i][j]*Qcon[j]; } @@ -1726,12 +1809,12 @@ u2p_solver_W(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose) 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[0]; - uint=pp[1]; + rho=pp[RHO]; + uint=pp[UU]; utcon[0]=0.; - utcon[1]=pp[2]; - utcon[2]=pp[3]; - utcon[3]=pp[4]; + utcon[1]=pp[VX]; + utcon[2]=pp[VY]; + utcon[3]=pp[VZ]; if (VELPRIM != VELR) { @@ -1740,9 +1823,6 @@ u2p_solver_W(ldouble *uu, ldouble *pp, void *ggg,int Etype,int verbose) ldouble qsq=0.; for(i=1;i<4;i++) -#ifdef APPLY_OMP_SIMD - //#pragma omp simd -#endif for(j=1;j<4;j++) qsq+=utcon[i]*utcon[j]*gg[i][j]; ldouble gamma2=1.+qsq; @@ -1990,6 +2070,197 @@ u2p_solver_Bonly(ldouble *uu, ldouble *pp, void *ggg) } // int u2p_solver_Bonly + +//********************************************************************** +//force-free solver +//********************************************************************** +int +u2p_solver_ff(ldouble *uu, ldouble *pp, void *ggg,int verbose) +{ + +#ifdef FORCEFREE + int i,j,k; + ldouble alpha,alphasq; + ldouble rho,entr,uint,gamma2,gamma,D,Sc,Bsq,Esq; + ldouble betavec[4]; + ldouble utcon[4],Econ[4],Ecov[4],Bcon[4],Bcov[4]; + ldouble Qcov[4],Qtildecov[4],Qtildecon[4]; + + /****************************/ + //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; + + ldouble pgamma=GAMMA; + #ifdef CONSISTENTGAMMA + pgamma=pick_gammagas(geom->ix,geom->iy,geom->iz); + #endif + + /****************************/ + //conserved quantities etc + + //alpha + alpha = geom->alpha; + alphasq = alpha*alpha; + + //beta vector components + betavec[0] = 0.; + betavec[1] = alphasq * GG[0][1]; + betavec[2] = alphasq * GG[0][2]; + betavec[3] = alphasq * GG[0][3]; + + //McKinney 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; + + indices_21(Bcon,Bcov,gg); + Bsq = dot(Bcon,Bcov); + + //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 + utcon[0] = 0.; + utcon[1] = Qtildecon[1]/Bsq; + utcon[2] = Qtildecon[2]/Bsq; + utcon[3] = Qtildecon[3]/Bsq; + + // get Lorentz factor + int ll,mm; + ldouble vsq=0.; + ldouble vsqmax = 1. - 1./(GAMMAMAXFF*GAMMAMAXFF); + int hitgammaceil=0; + for(ll=1;ll<4;ll++) + { + for(mm=1;mm<4;mm++) + { + vsq += gg[ll][mm]*utcon[ll]*utcon[mm]; + } + } + + ldouble vsq_orig = vsq; + + // limiter + if(vsq>vsqmax) + { + hitgammaceil=1; + ldouble vfac = sqrt(vsqmax/vsq); + utcon[1] *= vfac; + utcon[2] *= vfac; + utcon[3] *= vfac; + vsq = 0.; + for(ll=1;ll<4;ll++) + { + for(mm=1;mm<4;mm++) + { + vsq += gg[ll][mm]*utcon[ll]*utcon[mm]; + } + } + } + + gamma2 = 1./(1-vsq); + gamma = sqrt(gamma2); + + utcon[1] *= gamma; + utcon[2] *= gamma; + utcon[3] *= gamma; + + + // convert VELR to VELPRIM if necessary + if(VELPRIM != VELR) + { + conv_vels(utcon,utcon,VELR,VELPRIM,gg,GG); + } + + 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(utcon[1])) + { + if(verbose>0) printf("utcon[1] not finite in u2p_solver_ff %e %e %e \n",Bsq,Esq,gamma); getchar(); + return -104; + } + + + //************************************** + // Return new primitives + pp[VXFF]=utcon[1]; + pp[VYFF]=utcon[2]; + pp[VZFF]=utcon[3]; + + pp[B1]=uu[B1] * gdetu_inv; + pp[B2]=uu[B2] * gdetu_inv; + pp[B3]=uu[B3] * gdetu_inv; + + ////////////////////////////////////////////////////////////////////////////////////////// + // Hydro part + // ANDREW TODO -- replace with more sophisticated version + ///////////////////////////////////////////////////////////////////////////////////////// + + //************************************** + // Get density, entropy, and energy density + + D = uu[RHO] * alpha * gdetu_inv; // gdet rho ut, so D = gamma * rho + Sc = uu[ENTR] * alpha * gdetu_inv; // gdet S ut, so Sc = gamma * S + + rho = D/gamma; + entr = Sc/gamma; + 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; + } + + //*************************************** + //return new primitives + pp[RHO]=rho; + pp[UU]=uint; + pp[ENTR]=entr; + + // set velocities equal + pp[VX] = pp[VXFF]; + pp[VY] = pp[VYFF]; + pp[VZ] = pp[VZFF]; + + //if(verbose>0) printf("u2p_solver_ff returns 0\n"); + return 0; //ok +#endif +} + + //********************************************************************** //count the number of entropy inversions //********************************************************************** diff --git a/u2prad.c b/u2prad.c index 70a95e8..6ccef6d 100644 --- a/u2prad.c +++ b/u2prad.c @@ -553,7 +553,7 @@ check_floors_rad(ldouble *pp, int whichvel,void *ggg) //absolute rad-frame EE: if(pp[EE0]ix,geom->iy,geom->iz,pp[0],pp[EE0]); + if(verbose) printf("rhd_floors CASE R0 at (%d,%d,%d): %e %e\n",geom->ix,geom->iy,geom->iz,pp[RHO],pp[EE0]); pp[EE0]=ERADFLOOR; ret=-1; } @@ -561,32 +561,32 @@ check_floors_rad(ldouble *pp, int whichvel,void *ggg) #ifndef SKIPRADSOURCE //Ehat/rho ratios - if(Ehatix,geom->iy,geom->iz,pp[0],pp[EE0]); - pp[EE0]=Eratio*EERHORATIOMIN*pp[0]; + if(verbose) printf("hd_floors CASE R2 at (%d,%d,%d): %e %e\n",geom->ix,geom->iy,geom->iz,pp[RHO],pp[EE0]); + pp[EE0]=Eratio*EERHORATIOMIN*pp[RHO]; ret=-1; } - if(Ehat>EERHORATIOMAX*pp[0]) + if(Ehat>EERHORATIOMAX*pp[RHO]) { - if(verbose) printf("hd_floors CASE R3 at (%d,%d,%d): %e %e\n",geom->ix,geom->iy,geom->iz,pp[0],Ehat); - pp[0]=1./EERHORATIOMAX*Ehat; + if(verbose) printf("hd_floors CASE R3 at (%d,%d,%d): %e %e\n",geom->ix,geom->iy,geom->iz,pp[RHO],Ehat); + pp[RHO]=1./EERHORATIOMAX*Ehat; ret=-1; } //Ehat/uint ratios - if(Ehatix,geom->iy,geom->iz,pp[1],Ehat); - pp[EE0]=Eratio*EEUURATIOMIN*pp[1]; + if(verbose) printf("hd_floors CASE R4 at (%d,%d,%d): %e %e\n",geom->ix,geom->iy,geom->iz,pp[UU],Ehat); + pp[EE0]=Eratio*EEUURATIOMIN*pp[UU]; ret=-1; } - if(Ehat>EEUURATIOMAX*pp[1]) + if(Ehat>EEUURATIOMAX*pp[UU]) { - if(verbose) printf("hd_floors CASE R5 at (%d,%d,%d): %e %e\n",geom->ix,geom->iy,geom->iz,pp[1],Ehat); - pp[1]=1./EEUURATIOMAX*Ehat; + if(verbose) printf("hd_floors CASE R5 at (%d,%d,%d): %e %e\n",geom->ix,geom->iy,geom->iz,pp[UU],Ehat); + pp[UU]=1./EEUURATIOMAX*Ehat; ret=-1; }