diff --git a/PROBLEMS/RADSURVEY/bc.c b/PROBLEMS/RADSURVEY/bc.c new file mode 100755 index 0000000..295cf75 --- /dev/null +++ b/PROBLEMS/RADSURVEY/bc.c @@ -0,0 +1,288 @@ +//returns problem specific BC +//int calc_bc(int ix,int iy,int iz,ldouble t, +// ldouble *uu,ldouble *pp,int ifinit,int BCtype) + +/**********************/ +//geometries +ldouble gdet_src,gdet_bc,Fx,Fy,Fz,ppback[NV],ppback2[NV],pp1[NV],pp2[NV]; +int iix,iiy,iiz,iv; + +int iiytest=32; + +struct geometry geom; +fill_geometry(ix,iy,iz,&geom); + +/**********************/ + +#ifdef COPY_XBC //simple inflow/outflow +if(BCtype==XBCHI || BCtype==XBCLO) +{ + iiy=iy; + iiz=iz; + if(BCtype==XBCHI) iix=NX-1; + else if(BCtype==XBCLO) iix=0; + + for(iv=0;iv JET + +#define MINX 0 +#define MAXX 1. +#define Y_OFFSET 0.01 +#define MINY -(1.-Y_OFFSET) //10^-8 away from the poles seems to be the last safe point +#define MAXY 1.-Y_OFFSET + +#define MKSR0 0 //-1.35 // Should probably be 0 for jetcoords! (issue with ix=-2 metric) +#define HYPRBRK RMAX +#define FJET 0.4 +#define FDISK 0.4 + +#define RUNI RMIN +#define RCOLL_JET 10000 +#define RDECOLL_JET 2.//2*RMIN +#define RCOLL_DISK 10.//5*RMIN +#define RDECOLL_DISK 2.//2*RMIN + +#define ALPHA_1 1 +#define ALPHA_2 0.2//0.375 + +#define CYLINDRIFY +#define RCYL 20.//4.//10.//10 +#define NCYL 1.//0.5//1. +#endif + +#define PHIWEDGE (2*M_PI) +#define MINZ (-PHIWEDGE/2.) +#define MAXZ (PHIWEDGE/2.) + +/************************************/ +//boundary conditions +/************************************/ + +#define SPECIFIC_BC +#define PERIODIC_ZBC +//#define COPY_XBC //simpler outflowing boundary conditions in radius + +/************************************/ +//common physics +/************************************/ +#define GAMMA (5./3.) //(13./9.) +#define GAMMAI (5./3.) +#define GAMMAE (4./3.) + +#define HFRAC 1. //mass fraction of the hydrogen X +#define HEFRAC 0. //mass fraction of helium Y + +/************************************/ +//Initial Torus/Atmosphere +/************************************/ + +#define LIMOTORUS +#define NTORUS 3 //3 for SANE, 1 for MAD + +#ifdef LIMOTORUS +// Penna et all 2013 limotorus +#define LT_KAPPA 0.008 //This density normalization gives peak density ~1 + // in code units +#define LT_MAXDENS 1000 //This maximum density value will remove issues + // with high density points outside the disk + // we see in retrograde esp. + // be careful with scale,esp with LT_KAPPA! + +//LT_XI sets the torus size -- these values chosen to give outer radius ~500 +//a=0.75 0.752, a=0.5 0.757, a=0.25 0.761 +// a=-0.25 0.771, a=-.5 0.776, a=-.75 0.782 +#if(BHSPIN==-0.75) +#define LT_XI 0.782 +#elif(BHSPIN==-0.5) +#define LT_XI 0.776 +#elif(BHSPIN==-0.25) +#define LT_XI 0.771 +#if(BHSPIN==0.25) +#define LT_XI 0.761 +#elif(BHSPIN==0.5) +#define LT_XI 0.757 +#elif(BHSPIN==0.75) +#define LT_XI 0.752 +#endif + +#define LT_R1 42. //42. +#define LT_R2 800. //500. +#define LT_GAMMA 5./3. +#define LT_RIN 12. //10.5 + +#else + +// Fishbone-Moncrief torus +#define FM_rin 20. +#define FM_rmax 41. +#define FM_rho0 1.//rhoCGS2GU(1.e-18) //density normalization + +#endif + +// atmosphere //TODO better scalings! +#define ATMTYPE 0 +#define RHOATMMIN 1.e-4*pow(2.,-1.5) //*FM_rho0 +#define UINTATMMIN 1.e-6*pow(2.,-2.5)/(GAMMA-1.)/3. +#define ATMTRADINIT 2.7 +#define ERADATMMIN calc_LTE_EfromT(ATMTRADINIT) + +// B-field +#if(NTORUS==0) // single loop v1 +#define Aphi_rho_cut 0.2 +#undef MAXBETA +#define MAXBETA (100.) //target pgas/pmag inside torus +#define BETANORMFULL +#endif + +#if(NTORUS==1) // single loop v2 +#define Aphi_r_cut 400. +#define Aphi_rho_cut 0.2 +#undef MAXBETA +#define MAXBETA (100.) //target pgas/pmag inside torus +#define BETANORMFULL +#endif + +#if(NTORUS==2) // dipolar loops +#define Aphi_lambda 0.5 +#define Aphi_rchop 400. +#undef MAXBETA +#define MAXBETA (100.) //target pgas/pmag inside torus +#define BETANORMFULL +#endif + +#if(NTORUS==3) // quadrupolar loops +#define Aphi_lambda 1. +#define Aphi_rchop 400. +#undef MAXBETA +#define MAXBETA (100.) //target pgas/pmag inside torus +#define BETANORMFULL +#endif + +/************************************/ +//output +/************************************/ +//#define DTAVG .1 //how frequently to compute quantities included in avg +#define DTOUT1 10 //res +#define DTOUT2 100 //avg +#define DTOUT3 1000 //box,var + +//stopping condition + +#define NOUTSTOP 0 //1000//2000 + +//#define DUMPS_READ_HDF5 +//#define DUMPS_WRITE_HDF5 + +#define OUTCOORDS BLCOORDS //KSCOORDS +#define OUTVEL VEL4 + +#define COORDOUTPUT 0 +#define GRIDOUTPUT 0 +#define SILOOUTPUT 1 +#define OUTOUTPUT 0 +#define SIMOUTPUT 0 //needs to be nonzero for phiavg or phicorr +#define SIMOUTPUT_PHIAVG +//#define SIMOUTPUT_PHICORR +//#define CGSOUTPUT +//#define GRTRANSSIMOUTPUT_2 + +#define RADOUTPUT 0 +//#define RADOUTPUTWITHINDTHETA (M_PI/6.) + +#define SCAOUTPUT 0 +#define AVGOUTPUT 0 +#define NORELELAVGS +#define THOUTPUT 0 +//#define THPROFRADIUS 30 + + diff --git a/PROBLEMS/RADSURVEY/define_075.h b/PROBLEMS/RADSURVEY/define_075.h new file mode 100755 index 0000000..12aad9d --- /dev/null +++ b/PROBLEMS/RADSURVEY/define_075.h @@ -0,0 +1,421 @@ +/************************************/ +//general +/************************************/ +#define BHDISK_PROBLEMTYPE + +/************************************/ +//restart +/************************************/ +//#define NORESTART +#define RESTARTGENERALINDICES +#define RESTARTNUM -1 +//#define PERTURBUINT 0.02 + +/************************************/ +//blackhole +/************************************/ +#define MASS 6.5e9 +#define BHSPIN 0.75 +#define RH 1.66144 + +/************************************/ +//U2P-related choices +/************************************/ +//#define U2P_EQS U2P_EQS_JONS +//#define U2P_SOLVER U2P_SOLVER_WP +#define U2PCONV 1.e-10 + +/************************************/ +//magnetic choices +/************************************/ +#define MAGNFIELD +#define GDETIN 1 //must be 1 for MAGNFIELD +#define VECPOTGIVEN +#define INIT_MAGN_CORNERS //initialize magnetic field on corners/edges (which?) + +/************************************/ +//dynamo choices -- 2D ONLY +/************************************/ +/* +#define MIMICDYNAMO +#define CALCHRONTHEGO +#define THETAANGLE 0.25 +#define ALPHAFLIPSSIGN +#define ALPHADYNAMO 0.314 +#define DAMPBETA +#define BETASATURATED 0.1 +#define ALPHABETA 6.28 +*/ +/************************************/ +//radiation choices +/************************************/ +//#define RADIATION + +//#define SKIPRADSOURCE //advective only +//#define SKIPRADEVOLUTION //keeps initial values in place + +#ifdef RADIATION + +//#define BASICRADIMPLICIT +//#define RADIMPSTOPWHENFAIL +//#define RADIMP_START_WITH_BISECT +//#define BALANCEENTROPYWITHRADIATION + +//#define ALLOWRADCEILINGINIMPLICIT +//#define EVOLVEPHOTONNUMBER +#define OPDAMPINIMPLICIT 0 +#define SCALE_JACOBIAN +#endif + +//opacities +#define SCATTERING +#define BREMSSTRAHLUNG +//#define KLEINNISHINA +#define SYNCHROTRON +#define NO_SYNCHROTRON_BRIDGE_FUNCTIONS +#define MAXDIFFTRADS 1.e3 +#define MAXDIFFTRADSNEARBH 1.e2 + +//implicit convergence +#define RADIMPCONV 1.e-14 +#define RADIMPCONVREL 1.e-6 +#define RADIMPCONVRELERR (1.e-4) +#define RADIMPCONVRELENTR 1.e-6 +#define RADIMPCONVRELENTRERR 1.e-4 +#define RADIMPENTRCONV 1.e-5 +#define RADIMPEPS 1.e-8 +#define RADIMPMAXITER 100 +#define RADIMPLICITDAMPINGFACTOR 5. +#define RADIMPLICITMAXNPHCHANGE 100. +#define RADIMPLICITMAXENCHANGEDOWN 100. +#define RADIMPLICITMAXENCHANGEUP 10. +#define RADIMPLICITMAXTECHANGE 2. +#define IMPLICITMAXTGASCHANGE 2. +#define MAXRADIMPDAMPING 1.e-7 + +//Radiation Viscosity choices +#ifdef RADIATION + +#define RADVISCOSITY SHEARVISCOSITY +#define ACCELRADVISCOSITY +#define RADVISCMFPSPHMAX 10. +#define RADVISCMFPSPH +#define RADVISCNUDAMP +#define RADVISCMAXVELDAMP +#define ALPHARADVISC 0.1 +#define MAXRADVISCVEL 0.1 + +//#define SKIPHDEVOLUTION +//#define SKIPRADEVOLUTION +//#define SKIPEVOLUTION +//#define SKIPRADSOURCE +//#define SKIPCOULOMBCOUPLING +//#define RADIMPSTOPWHENFAIL +#define DAMPCOMPTONIZATIONATBH +//#define NO_COMPTONIZATION +//#define SKIPFANCYOPACITIES +//#define ENFORCEENTROPY +//#define GASRADCOUPLEDWAVESPEEDS + +#endif + +/************************************/ +//electron choices +/************************************/ +//#define EVOLVEELECTRONS + +#ifdef EVOLVEELECTRONS +#define CONSISTENTGAMMA +#define GAMMAINTCONSISTENTWITHCV //Ramesh's routine for inverting gamma_int +#define RELELENTROPY + +//heating +#define HEATELECTRONS +//#define HEATELECTRONS_HOWES +//#define HEATELECTRONS_ROWAN +#define HEATELECTRONS_ROWAN2 +//#define HEATELECTRONS_ROWAN3 + +#define NOHEATATBH + +//#define HEATELECTRONSATENDRK2 +//#define DISSIPATIONFROMGASONLY +#define FORCEGAMMAGASFIXED + +//entropy mixing +//#define MIXENTROPIESPROPERLY +//#define UPWINDENTROPYMIXING +//#define DONOTLIMITENTRINMIXING + +//silo output +//#define PRINTVISCHEATINGTOSILO +//#define PRINTCOULOMBTOSILO + +//floors +#define UEUINTMINRATIO 1.e-3 +#define UIUINTMINRATIO 1.e-3 +#define TEMPEMINIMAL 1.e2 +#define TEMPIMINIMAL 1.e2 +#define TEMPEMINIMALFRACTION 1.e-6 +#define TEMPIMINIMALFRACTION 1.e-6 +#define TEMPEMAXIMALFRACTION 1.e2 +#define TEMPIMAXIMALFRACTION 1.e2 +#endif + +/************************************/ +//reconstruction / Courant +/************************************/ +#define INT_ORDER 2 +#define TIMESTEPPING RK2IMEX +#define TSTEPLIM .9 +#define FLUXMETHOD LAXF_FLUX + +#define FLUXLIMITER 0 +#define MINMOD_THETA 1.9 +#define SHUFFLELOOPS 0 + +#define DOFIXUPS 1 +#define DOU2PRADFIXUPS 0 //ANDREW -- these fixups lead to failures! +#define DOU2PMHDFIXUPS 1 +#define DORADIMPFIXUPS 1 + +/************************************/ +//rmhd floors +/************************************/ +#define CORRECT_POLARAXIS +#define NCCORRECTPOLAR 2 +#define UURHORATIOMIN 1.e-8 +#define UURHORATIOMAX 1.e2 +#define EERHORATIOMIN 1.e-20 +#define EERHORATIOMAX 1.e20 +#define EEUURATIOMIN 1.e-20 +#define EEUURATIOMAX 1.e20 +#define B2UURATIOMIN 0. +#define B2UURATIOMAX 1.e3 +#define B2RHORATIOMIN 0. +#define B2RHORATIOMAX 100. +#define GAMMAMAXRAD 20. +#define GAMMAMAXHD 20. +#define RHOFLOOR 1.e-40 + +/************************************/ +//resolution +/************************************/ +//total resolution +#define TNX 320 //336//256//32//128//312 +#define TNY 256 //192//192//32//128//200 +#define TNZ 1 //128//192//192 + +#define SILO2D_XZPLANE + +//number of tiles +#define NTX 20 //28 +#define NTY 16 //16 +#define NTZ 1//8 //16 + +/************************************/ +//coordinates +/************************************/ +//#define myJETCOORDS +#define myMKS2COORDS +//#define myMKS3COORDS + +#define METRICAXISYMMETRIC + +#define RMIN 0.825*RH +#define RMAX 1.e5 + +#ifdef myMKS2COORDS //modified Kerr-Shild +#define METRICNUMERIC +#define MYCOORDS MKS2COORDS + +#define MKSR0 0.//35 +#define MKSH0 0.7 +#define MKSMY1 0.002 +#define MKSMY2 0.02 +#define MKSMP0 1.3 + +#define MINX (log(RMIN-MKSR0)) +#define MAXX (log(RMAX-MKSR0)) +#define MINY (0.01) +#define MAXY (1. - 0.01) + +#endif + +#ifdef myMKS3COORDS //modified Kerr-Shild further from axis +#define METRICNUMERIC +#define MYCOORDS MKS3COORDS +#define MINX (log(RMIN-MKSR0)) +#define MAXX (log(RMAX-MKSR0)) +#define MINY 0. +#define MAXY 1. + +#define MKSR0 -1.35 +#define MKSH0 0.7 +#define MKSMY1 0.002 +#define MKSMY2 0.02 +#define MKSMP0 1.3 +#endif + +#ifdef myJETCOORDS //concentrate resolution in jet and disk zones +#define MYCOORDS JETCOORDS + +#define METRICNUMERIC +//#define DERIVS_NOGSL // use a faster numeric derivative for coordinate transformations +#define PRECOMPUTE_MY2OUT // precompute transformation matrices from BL <--> JET + +#define MINX 0 +#define MAXX 1. +#define Y_OFFSET 0.001 +#define MINY -(1.-Y_OFFSET) //10^-8 away from the poles seems to be the last safe point +#define MAXY 1.-Y_OFFSET + +#define MKSR0 0 //-1.35 // Should probably be 0 for jetcoords! (issue with ix=-2 metric) +#define HYPRBRK 1000 +#define FJET 0.25 +#define FDISK 0.4 + +//#define RUNI RMIN +//#define RCOLL_JET 1000 +//#define RDECOLL_JET 2*RMIN//2*RMIN +//#define RCOLL_DISK 20*RMIN//5*RMIN +//#define RDECOLL_DISK 2*RMIN//2*RMIN + +#define RUNI RMIN +#define RCOLL_JET 1000 +#define RDECOLL_JET 2*RMIN +#define RCOLL_DISK 5*RMIN +#define RDECOLL_DISK 2*RMIN + +#define ALPHA_1 1 +#define ALPHA_2 0.375 + +#define CYLINDRIFY +#define RCYL 40.//4.//10.//10 +#define NCYL 1.//0.5//1. +#endif + +#define PHIWEDGE (2*M_PI) +#define MINZ (-PHIWEDGE/2.) +#define MAXZ (PHIWEDGE/2.) + +/************************************/ +//boundary conditions +/************************************/ + +#define SPECIFIC_BC +#define PERIODIC_ZBC +//#define COPY_XBC //simpler outflowing boundary conditions in radius + +/************************************/ +//common physics +/************************************/ +#define GAMMA (5./3.) //(13./9.) +#define GAMMAI (5./3.) +#define GAMMAE (4./3.) + +#define HFRAC 1. //mass fraction of the hydrogen X +#define HEFRAC 0. //mass fraction of helium Y + +/************************************/ +//Initial Torus/Atmosphere +/************************************/ + +#define LIMOTORUS +#define NTORUS 1 + +#ifdef LIMOTORUS +// Penna et all 2013 limotorus +#define LT_KAPPA 0.008 //3.e8 //TODO what should the density normalization be? +#define LT_XI 0.752 //0.7545 //0.7135 +#define LT_R1 42. //42. +#define LT_R2 800. //500. +#define LT_GAMMA 5./3. +#define LT_RIN 12. //10.5 + +#else + +// Fishbone-Moncrief torus +#define FM_rin 20. +#define FM_rmax 41. +#define FM_rho0 1.//rhoCGS2GU(1.e-18) //density normalization + +#endif + +// atmosphere //TODO better scalings! +#define ATMTYPE 0 +#define RHOATMMIN 1.e-4*pow(2.,-1.5) //*FM_rho0 +#define UINTATMMIN 1.e-6*pow(2.,-2.5)/(GAMMA-1.)/3. +#define ATMTRADINIT 2.7 +#define ERADATMMIN calc_LTE_EfromT(ATMTRADINIT) + +// B-field +#if(NTORUS==0) // single loop v1 +#define Aphi_rho_cut 0.2 +#undef MAXBETA +#define MAXBETA (100.) //target pgas/pmag inside torus +#define BETANORMFULL +#endif + +#if(NTORUS==1) // single loop v2 +#define Aphi_r_cut 400. +#define Aphi_rho_cut 0.2 +#undef MAXBETA +#define MAXBETA (100.) //target pgas/pmag inside torus +#define BETANORMFULL +#endif + +#if(NTORUS==2) // dipolar loops +#define Aphi_lambda 0.5 +#define Aphi_rchop 200. +#undef MAXBETA +#define MAXBETA (100.) //target pgas/pmag inside torus +#define BETANORMFULL +#endif + +#if(NTORUS==3) // quadrupolar loops +#define Aphi_lambda 1. +#define Aphi_rchop 200. +#undef MAXBETA +#define MAXBETA (100.) //target pgas/pmag inside torus +#define BETANORMFULL +#endif + +/************************************/ +//output +/************************************/ +//#define DTAVG .1 //how frequently to compute quantities included in avg +#define DTOUT1 10 //res +#define DTOUT2 100 //avg +#define DTOUT3 1000 //box,var + +//stopping condition + +#define NOUTSTOP 0 //1000//2000 + +//#define DUMPS_READ_HDF5 +//#define DUMPS_WRITE_HDF5 + +#define OUTCOORDS KSCOORDS//BLCOORDS +#define OUTVEL VEL4 + +#define COORDOUTPUT 0 +#define GRIDOUTPUT 0 +#define SILOOUTPUT 1 +#define OUTOUTPUT 0 +#define SIMOUTPUT 0 //needs to be nonzero for phiavg or phicorr +#define SIMOUTPUT_PHIAVG +//#define SIMOUTPUT_PHICORR +//#define CGSOUTPUT +//#define GRTRANSSIMOUTPUT_2 + +#define RADOUTPUT 0 +//#define RADOUTPUTWITHINDTHETA (M_PI/6.) + +#define SCAOUTPUT 0 +#define AVGOUTPUT 0 +#define NORELELAVGS +#define THOUTPUT 0 +//#define THPROFRADIUS 30 + + diff --git a/PROBLEMS/RADSURVEY/init.c b/PROBLEMS/RADSURVEY/init.c new file mode 100644 index 0000000..7afebe4 --- /dev/null +++ b/PROBLEMS/RADSURVEY/init.c @@ -0,0 +1,340 @@ +int init_dsandvels_limotorus(FTYPE r, FTYPE th, FTYPE a, FTYPE *rhoout, FTYPE *uuout, FTYPE *ell); +ldouble rhomax_limotorus_approx(FTYPE a); +int init_dsandvels_fishbone_moncrief(FTYPE r, FTYPE th, FTYPE a, FTYPE *rhoout, FTYPE *uuout, FTYPE *ell); + +ldouble rho,mx,my,mz,m,E,uint,pgas,Fx,Fy,Fz,pLTE,ell; +ldouble uu[NV], pp[NV],ppback[NV],T,uintorg; +ldouble Vphi,Vr; +ldouble D,W,eps,uT,uphi,uPhi; + +//geometries +struct geometry geom; +fill_geometry(ix,iy,iz,&geom); + +struct geometry geomBL; +fill_geometry_arb(ix,iy,iz,&geomBL,KERRCOORDS); + +ldouble r=geomBL.xx; +ldouble th=geomBL.yy; + +// good faith estimate of initial gamma_gas +#ifdef CONSISTENTGAMMA +ldouble gamma=calc_gammaintfromTei(1.e10,1.e10); +set_u_scalar(gammagas,ix,iy,iz,gamma); +#endif + +// calculate torus rho, uint, ell from tools.c function +#ifdef LIMOTORUS +init_dsandvels_limotorus(r, th, BHSPIN, &rho, &uint, &ell); //penna+ torus +//printf("%d %d %d %.2e\n",ix,iy,iz,rho); +#else +init_dsandvels_fishbone_moncrief(r, th, BHSPIN, &rho, &uint, &ell); //fishbone-moncreif +#endif +uintorg=uint; + +// maximum density -- removes issues with cells outside the disk +#ifdef LT_MAXDENS +if(rho>LT_MAXDENS) + rho=-1.; +#endif + +if(rho<0.) // we are outside the donut, in the atmosphere +{ + // atmosphere: atmtype=0. + // Need to define RHOATMMIN, UINTATMMIN defined at rout=2. + set_hdatmosphere(pp,geom.xxvec,geom.gg,geom.GG,ATMTYPE); //density, temperature + #ifdef RADIATION + set_radatmosphere(pp,geom.xxvec,geom.gg,geom.GG,0); //radiation + #ifdef EVOLVEPHOTONNUMBER + pp[NF]=1./2.70118/K_BOLTZ * pp[EE] / ATMTRADINIT; //atmosphere has constant init temp. + #endif + #endif +} +else // we are inside the donut +{ + + //calculate the hydro atmosphere values as backup + set_hdatmosphere(ppback,geom.xxvec,geom.gg,geom.GG,ATMTYPE); + #ifdef RADIATION + //the atmosphere radiation remains as background + set_radatmosphere(pp,geom.xxvec,geom.gg,geom.GG,0); + #endif + + if(ppback[0]>rho) + printf("INSIDE LARGE ATM: %.4e %.4e \n", rho, ppback[0]); + + //Calculate torus 4-velocity from ell; + ldouble ult,ulph,ucov[4],ucon[4]; + + // assuming ell = u_phi / |u_t| + #ifdef LIMOTORUS + + ell*=-1.; + ulph = sqrt(-1./(geomBL.GG[0][0]/ell/ell + 2./ell*geomBL.GG[0][3] + geomBL.GG[3][3])); + ult = ulph / ell; + + ucov[0]=ult; + ucov[1]=0.; + ucov[2]=0.; + ucov[3]=ulph; + indices_12(ucov,ucon,geomBL.GG); + + #else + + // Calculate torus 4-velocity from ell; new method + // Based on the statement in Kozlowski et al. that ell is actually u_phi u^t + ldouble rhosq = geomBL.gg[0][0] * geomBL.gg[3][3] - geomBL.gg[0][3] * geomBL.gg[0][3]; + ldouble utsq = (-geomBL.gg[3][3] - sqrt(geomBL.gg[3][3] * geomBL.gg[3][3] - 4. * rhosq * ell * ell)) / (2. * rhosq); + + ucon[0] = sqrt(utsq); + ucon[1] = 0.; + ucon[2] = 0.; + ucon[3] = ell / (geomBL.gg[3][3] * ucon[0]) - geomBL.gg[0][3] * ucon[0] / geomBL.gg[3][3]; + indices_21(ucon,ucov,geomBL.gg); + + #endif //LIMOTORUS + + conv_vels_ut(ucon,ucon,VEL4,VELPRIM,geomBL.gg,geomBL.GG); + + // set pp[0], pp[1], but use the atmosphere values as floors + //uint = LT_KAPPA * pow(rho, LT_GAMMA) / (LT_GAMMA - 1.); + //pgas = GAMMAM1 * uint; + pp[0]=my_max(rho,ppback[0]); //choice of higher of two densities between background and torus + pp[1]=my_max(uint,ppback[1]); + + // set initial velocities + pp[2]=ucon[1]; + pp[3]=ucon[2]; + pp[4]=ucon[3]; + +#ifdef MAGNFIELD + pp[B1]=0.; + pp[B2]=0.; + pp[B3]=0.; //initially set B zero not to break coordinate transformation +#endif + +#ifdef RADIATION + pp[FX]=pp[VX]; //make the initial radiation atm in torus rotate with the same velocity + pp[FY]=pp[VY]; + pp[FZ]=pp[VZ]; +#ifdef EVOLVEPHOTONNUMBER + pp[NF]=1./2.70118/K_BOLTZ * pp[EE]/ATMTRADINIT; //initial rad atm has constant temp +#endif +#endif + + //transforming primitives from BL to MYCOORDS +#ifdef PRECOMPUTE_MY2OUT + trans_pall_coco_out2my(pp,pp,&geomBL,&geom); // Take advantage of precomputed MY <--> BL +#else + trans_pall_coco(pp, pp, BLCOORDS, MYCOORDS, geomBL.xxvec,&geomBL,&geom); +#endif + + //Set vector potential to compute magnetic field +#ifdef MAGNFIELD + + ldouble Acov[4]; + ldouble r_mag,th_mag; + Acov[0]=Acov[1]=Acov[2]=0.; + +#ifdef INIT_MAGN_CORNERS + + // ANDREW define at corners not cell centers + // ANDREW -- shouldn't this actually be edges? + // TODO: What about boundary corners? + // TODO: We assume we are well inside grid so that A=0 out there! + 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 + +#if (NTORUS==0) // single loop v1 + + //recompute rho (possibly at corner) and torus max rho + ldouble rhomax; + #ifdef LIMOTORUS + init_dsandvels_limotorus(r_mag, th_mag, BHSPIN, &rho, &uint, &ell); + rhomax = rhomax_limotorus_approx(BHSPIN); + #else + init_dsandvels_fishbone_moncrief(r_mag, th_mag, BHSPIN, &rho, &uint, &ell); + rhomax = FM_rho0; + #endif + if(fabs((rho-pp[RHO])/pp[RHO])>.5) + rho=pp[RHO]; + + //standard single poloidal loop: Aphi = max[(rho/rhomax) - Aphi_rho_cut, 0] + Acov[3]=my_max((rho/rhomax) - Aphi_rho_cut, 0.); + +#elif (NTORUS==1) // single loop v2 -- used for MAD code comparison + //standard single poloidal loop: Aphi = max[(rho/rhomax) - FM_Aphi_cut, 0] + + //recompute rho (possibly at corner) and torus max rho + ldouble rhomax,rin; + #ifdef LIMOTORUS + init_dsandvels_limotorus(r_mag, th_mag, BHSPIN, &rho, &uint, &ell); + rhomax = rhomax_limotorus_approx(BHSPIN); + rin = LT_RIN; + #else + init_dsandvels_fishbone_moncrief(r_mag, th_mag, BHSPIN, &rho, &uint, &ell); + rhomax = FM_rho0; + rin = FM_rin; + #endif + if(fabs((rho-pp[RHO])/pp[RHO])>.5) + rho=pp[RHO]; + + //vector potential + ldouble q; + q = (rho/rhomax)*pow(r_mag/rin, 3)*pow(sin(th_mag),3)*exp(-r_mag / Aphi_r_cut) - Aphi_rho_cut; + Acov[3]=my_max(q, 0.); + +#elif (NTORUS==2) // dipolar loops + ldouble lambda = Aphi_lambda; + ldouble anorm = 1.; + ldouble rchop = Aphi_rchop; //outer boundary of field loops + ldouble u_av, u_av_chop, u_av_mid; + ldouble rin; + + #ifdef LIMOTORUS + init_dsandvels_limotorus(r_mag, th_mag, BHSPIN, &rho, &u_av, &ell); // current loc + init_dsandvels_limotorus(r_mag, M_PI/2., BHSPIN, &rho, &u_av_mid, &ell); // midplane + init_dsandvels_limotorus(rchop, M_PI/2., BHSPIN, &rho, &u_av_mid, &ell); // midplane at rchop + rin = LT_RIN; + #else + init_dsandvels_fishbone_moncrief(r_mag, th_mag, BHSPIN, &rho, &u_av, &ell); // current loc + init_dsandvels_fishbone_moncrief(r_mag, M_PI/2., BHSPIN, &rho, &u_av_mid, &ell); // midplane + init_dsandvels_fishbone_moncrief(rchop, M_PI/2., BHSPIN, &rho, &u_av_mid, &ell); // midplane at rchop + rin = FM_rin; + #endif + + //vector potential follows contours of UU + ldouble uchop = u_av - u_av_chop; //vpot->zero on contour of radius r=rchop + ldouble uchopmid = u_av_mid - u_av_chop; //vpot->zero away from midplane + + ldouble STARTFIELD = 1.25*rin; + ldouble q, fr, fr_start, vpot=0.; + if (r_mag > STARTFIELD && r_mag < rchop) + { + q = anorm * (uchop - 0.2*uchopmid) / (0.8*uchopmid) * pow(sin(th_mag), 3); + } + else + q = 0; + + if(q > 0.) + { + fr = (pow(r_mag,0.6)/0.6 + 0.5/0.4*pow(r_mag,-0.4)) / lambda; + fr_start = (pow(STARTFIELD,0.6)/0.6 + 0.5/0.4*pow(STARTFIELD,-0.4)) / lambda; + vpot += q * sin(fr - fr_start) ; + } + + Acov[3]=vpot; + +#elif (NTORUS==3) // quadrupolar loops + ldouble lambda = Aphi_lambda; + ldouble anorm = 1.; + ldouble rchop = Aphi_rchop; //outer boundary of field loops + ldouble u_av, u_av_chop, u_av_mid; + ldouble rin; + + #ifdef LIMOTORUS + init_dsandvels_limotorus(r_mag, th_mag, BHSPIN, &rho, &u_av, &ell); // current loc + init_dsandvels_limotorus(r_mag, M_PI/2., BHSPIN, &rho, &u_av_mid, &ell); // midplane + init_dsandvels_limotorus(rchop, M_PI/2., BHSPIN, &rho, &u_av_mid, &ell); // midplane at rchop + rin = LT_RIN; + #else + init_dsandvels_fishbone_moncrief(r_mag, th_mag, BHSPIN, &rho, &u_av, &ell); // current loc + init_dsandvels_fishbone_moncrief(r_mag, M_PI/2., BHSPIN, &rho, &u_av_mid, &ell); // midplane + init_dsandvels_fishbone_moncrief(rchop, M_PI/2., BHSPIN, &rho, &u_av_mid, &ell); // midplane at rchop + rin = FM_rin; + #endif + + //vector potential follows contours of UU + ldouble uchop = u_av - u_av_chop; //vpot->zero on contour of radius r=rchop + ldouble uchopmid = u_av_mid - u_av_chop; //vpot->zero away from midplane + + ldouble STARTFIELD = 1.25*rin; + ldouble q, fr, fr_start, vpot=0.; + if (r_mag > STARTFIELD && r_mag < rchop) + { + q = anorm * (uchop - 0.2*uchopmid) / (0.8*uchopmid) * pow(sin(th_mag), 3); + } + else + q = 0; + + if(q > 0.) + { + fr = (pow(r_mag,0.6)/0.6 + 0.5/0.4*pow(r_mag,-0.4)) / lambda; + fr_start = (pow(STARTFIELD,0.6)/0.6 + 0.5/0.4*pow(STARTFIELD,-0.4)) / lambda; + vpot += q * sin(fr - fr_start) ; + } + + vpot *= sin((M_PI/2.-th_mag)); // quadrupolar part + + Acov[3]=vpot; +#endif + + // 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 + +} //if rho<0 + +// Entropy +pp[ENTR]=calc_Sfromu(pp[RHO],pp[UU],ix,iy,iz); + +// Thermodynamics +#ifdef EVOLVEELECTRONS + +ldouble rhogas=pp[RHO]; +ldouble Tgas=calc_PEQ_Tfromurho(pp[UU],pp[RHO],ix,iy,iz); + +//slightly colder electrons initially +ldouble ue=1./20.*pp[UU]; +ldouble ui=(1.-1./20.)*pp[UU]; +ldouble Te,Ti; +pp[ENTRE]=calc_Sefromrhou(calc_thermal_ne(pp)*MU_E*M_PROTON,ue,ELECTRONS); +pp[ENTRI]=calc_Sefromrhou(rhogas,ui,IONS); + +#ifdef CONSISTENTGAMMA +Ti=solve_Teifromnmu(pp[RHO]/MU_I/M_PROTON, M_PROTON, ui, IONS); +Te=solve_Teifromnmu(pp[RHO]/MU_E/M_PROTON, M_ELECTR, ue, ELECTRONS); +gamma=calc_gammaintfromTei(Te,Ti); //recompute gamma_gas from actual electron/ion temps +set_u_scalar(gammagas,ix,iy,iz,gamma); +#endif + +pp[ENTRE]=calc_SefromrhoT(rhogas,Te,ELECTRONS); +pp[ENTRI]=calc_SefromrhoT(rhogas,Ti,IONS); + +#endif //EVOLVEELECTRONS + +//all primitives to conserved +p2u(pp,uu,&geom); + +/***********************************************/ +// set the final values +/***********************************************/ + +int iv; + +for(iv=0;iv BL +#ifdef PRECOMPUTE_MY2OUT + trans_pall_coco_out2my(pp,pp,&geomBL,&geom); +#else + trans_pall_coco(pp, pp, BLCOORDS, MYCOORDS, geomBL.xxvec,&geomBL,&geom); +#endif + + //Set vector potential to compute magnetic field +#ifdef MAGNFIELD + + ldouble Acov[4]; + ldouble r_mag,th_mag; + Acov[0]=Acov[1]=Acov[2]=0.; + +#ifdef INIT_MAGN_CORNERS + + // ANDREW define at corners not cell centers + // ANDREW -- shouldn't this actually be edges? + // What about boundary corners? -- we assume we are well inside grid so that A=0 out there! + + 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 + +#if (NTORUS==0) // single loop v1 + //standard single poloidal loop: Aphi = max[(rho/rhomax) - FM_Aphi_cut, 0] + //recompute rho -- possibly at corner + + init_dsandvels_fishbone_moncrief(r_mag, th_mag, BHSPIN, &rho, &uint, &ell); + if(fabs((rho-pp[RHO])/pp[RHO])>.5) + rho=pp[RHO]; + + //vector potential + Acov[3]=my_max((rho/FM_rho0) - FM_Aphi_cut, 0.); + +#elif (NTORUS==1) // single loop v2 -- MAD code comparison + //standard single poloidal loop: Aphi = max[(rho/rhomax) - FM_Aphi_cut, 0] + //recompute rho -- possibly at corner + + init_dsandvels_fishbone_moncrief(r_mag, th_mag, BHSPIN, &rho, &uint, &ell); + if(fabs((rho-pp[RHO])/pp[RHO])>.5) + rho=pp[RHO]; + + ldouble q = (rho/FM_rho0)*pow(r_mag/FM_rin, 3); + q *= pow(sin(th_mag),3); + q *= exp(-r_mag / FM_rcut); + q -= FM_Aphi_cut; + + //vector potential + Acov[3]=my_max(q, 0.); + +#elif (NTORUS==2) // dipolar loops + ldouble lambda = 0.5; + ldouble anorm = 1.; + ldouble rchop = 90.; + ldouble u_av = uintorg; + ldouble u_av_chop, u_av_mid; + + //midplane at r_mag + init_dsandvels_fishbone_moncrief(r_mag, M_PI/2., BHSPIN, &rho, &u_av_mid, &ell); + //midplane at rchop + init_dsandvels_fishbone_moncrief(rchop, M_PI/2., BHSPIN, &rho, &u_av_mid, &ell); + //vector potential follows contours of UU + ldouble uchop = u_av - u_av_chop; //vpot->zero on contour of radius r=rchop + ldouble uchopmid = u_av_mid - u_av_chop; //vpot->zero away from midplane + + + ldouble rin=FM_rin; + ldouble STARTFIELD = 1.25*rin; + ldouble q, fr, fr_start, vpot=0.; + if (r_mag > STARTFIELD && r_mag < rchop) + { + q = anorm * (uchop - 0.2*uchopmid) / (0.8*uchopmid) * pow(sin(th_mag), 3); + } + else + q = 0; + + if(q > 0.) + { + fr = (pow(r_mag,0.6)/0.6 + 0.5/0.4*pow(r_mag,-0.4)) / lambda; + fr_start = (pow(STARTFIELD,0.6)/0.6 + 0.5/0.4*pow(STARTFIELD,-0.4)) / lambda; + vpot += q * sin(fr - fr_start) ; + } + + Acov[3]=vpot + +#elif (NTORUS==3) // quadrupolar loops + ldouble lambda = 0.5; + ldouble anorm = 1.; + ldouble rchop = 90.; + ldouble u_av = uintorg; + ldouble u_av_chop, u_av_mid; + + //midplane at r_mag + init_dsandvels_fishbone_moncrief(r_mag, M_PI/2., BHSPIN, &rho, &u_av_mid, &ell); + //midplane at rchop + init_dsandvels_fishbone_moncrief(rchop, M_PI/2., BHSPIN, &rho, &u_av_mid, &ell); + //vector potential follows contours of UU + ldouble uchop = u_av - u_av_chop; //vpot->zero on contour of radius r=rchop + ldouble uchopmid = u_av_mid - u_av_chop; //vpot->zero away from midplane + + ldouble rin=FM_rin; + ldouble STARTFIELD = 1.25*rin; + ldouble q, fr, fr_start, vpot=0.; + if (r_mag > STARTFIELD && r_mag < rchop) + { + q = anorm * (uchop - 0.2*uchopmid) / (0.8*uchopmid) * pow(sin(th_mag), 3); + } + else + q = 0; + + if(q > 0.) + { + fr = (pow(r_mag,0.6)/0.6 + 0.5/0.4*pow(r_mag,-0.4)) / lambda; + fr_start = (pow(STARTFIELD,0.6)/0.6 + 0.5/0.4*pow(STARTFIELD,-0.4)) / lambda; + vpot += q * sin(fr - fr_start) ; + } + + Acov[3]=vpot*sin((M_PI/2.-th_mag));; +#endif + + // 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 + +} //if rho<0 + +// Entropy +pp[ENTR]=calc_Sfromu(pp[RHO],pp[UU],ix,iy,iz); + +// Thermodynamics +#ifdef EVOLVEELECTRONS + +ldouble rhogas=pp[RHO]; +ldouble Tgas=calc_PEQ_Tfromurho(pp[UU],pp[RHO],ix,iy,iz); + +//slightly colder electrons initially +ldouble ue=1./20.*pp[UU]; +ldouble ui=(1.-1./20.)*pp[UU]; +ldouble Te,Ti; +pp[ENTRE]=calc_Sefromrhou(calc_thermal_ne(pp)*MU_E*M_PROTON,ue,ELECTRONS); +pp[ENTRI]=calc_Sefromrhou(rhogas,ui,IONS); + +#ifdef CONSISTENTGAMMA +Ti=solve_Teifromnmu(pp[RHO]/MU_I/M_PROTON, M_PROTON, ui, IONS); +Te=solve_Teifromnmu(pp[RHO]/MU_E/M_PROTON, M_ELECTR, ue, ELECTRONS); + + +gamma=calc_gammaintfromTei(Te,Ti); //recompute gamma_gas from actual electron/ion temps +set_u_scalar(gammagas,ix,iy,iz,gamma); +#endif + +pp[ENTRE]=calc_SefromrhoT(rhogas,Te,ELECTRONS); +pp[ENTRI]=calc_SefromrhoT(rhogas,Ti,IONS); + +#endif //EVOLVEELECTRONS + +//all primitives to conserved +p2u(pp,uu,&geom); + +/***********************************************/ +// set the final values +/***********************************************/ + +int iv; + +for(iv=0;iv BL + #ifdef PRECOMPUTE_MY2OUT + trans_pall_coco_out2my(pp,pp,&geomBL,&geom); + #else + trans_pall_coco(pp, pp, BLCOORDS, OUTCOORDS, geomBL.xxvec,&geomBL,&geom); + #endif + + //Set Magnetic Field from vector potential +#ifdef MAGNFIELD + //MYCOORDS vector potential to calculate B's + ldouble Acov[4]; + ldouble r_mag,th_mag; + Acov[0]=Acov[1]=Acov[2]=0.; + +#ifdef INIT_MAGN_CORNERS + + //ANDREW define at corners not cell centers + //ANDREW -- shouldn't this actually be edges? + //what about boundary corners? -- assume we are well inside grid so that A=0 out there? + 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]; + + //define bfield at corners not cell centers + r_mag=r_c; + th_mag=th_c; +#else + r_mag=r; + th_mag=th; +#endif + +#if (NTORUS==4) + //LIMOFIELD from a=0 SANE harm init.c + denser dipolar loops + ldouble lambda = 1.75; //ANDREW changed from 1 + ldouble anorm=1.; //BOBMARK: not used, letting HARM normalize the field + ldouble rchop = 550.; //outer boundary of field loops + ldouble u_av = pp[UU]; + ldouble u_av_chop, u_av_mid; + //midplane at r_mag + init_dsandvels_limotorus(r_mag, M_PI/2., BHSPIN, &rho, &u_av_mid, &ell); + //midplane at rchop + init_dsandvels_limotorus(rchop, M_PI/2., BHSPIN, &rho, &u_av_chop, &ell); + + //vetor potential follows contours of UU + ldouble uchop = u_av - u_av_chop; //vpot->zero on contour of radius r=rchop + ldouble uchopmid = u_av_mid - u_av_chop; //vpot->zero away from midplane + + ldouble rin=LT_RIN; + ldouble STARTFIELD = 2.*rin; + ldouble q, fr, fr_start, vpot=0.; + if (r_mag > STARTFIELD && r_mag < rchop) { + q = anorm * (uchop - 0.2*uchopmid) / (0.8*uchopmid) * pow(sin(th_mag), 3); // * pow(tanh(r/rsmooth),2); + } else q = 0; + + if(q > 0.) { + fr = (pow(r_mag,0.6)/0.6 + 0.5/0.4*pow(r_mag,-0.4)) / lambda; + fr_start = (pow(STARTFIELD,0.6)/0.6 + 0.5/0.4*pow(STARTFIELD,-0.4)) / lambda; + vpot += q * sin(fr - fr_start) ; + } + + Acov[3] = vpot; + +#elif (NTORUS==10) + //LIMOFIELD from a=0 SANE harm init.c + denser quadrupolar loops + ldouble lambda = 1.; + ldouble anorm=1.; //BOBMARK: not used, letting HARM normalize the field + ldouble rchop = 550.; //outer boundary of field loops + ldouble u_av = pp[UU]; + ldouble u_av_chop, u_av_mid; + //midplane at r_mag + init_dsandvels_limotorus(r_mag, M_PI/2., BHSPIN, &rho, &u_av_mid, &ell); + //midplane at rchop + init_dsandvels_limotorus(rchop, M_PI/2., BHSPIN, &rho, &u_av_chop, &ell); + + //vector potential follows contours of UU + ldouble uchop = u_av - u_av_chop; //vpot->zero on contour of radius r=rchop + ldouble uchopmid = u_av_mid - u_av_chop; //vpot->zero away from midplane + + ldouble rin=LT_RIN; + ldouble STARTFIELD = 2.*rin; + ldouble q, fr, fr_start, vpot=0.; + if (r_mag > STARTFIELD && r_mag < rchop) { + q = anorm * (uchop - 0.2*uchopmid) / (0.8*uchopmid) * pow(sin(th_mag), 3); + } else q = 0; + + if(q > 0.) { + fr = (pow(r_mag,0.6)/0.6 + 0.5/0.4*pow(r_mag,-0.4)) / lambda; + fr_start = (pow(STARTFIELD,0.6)/0.6 + 0.5/0.4*pow(STARTFIELD,-0.4)) / lambda; + vpot += q * sin(fr - fr_start) ; + } + + Acov[3] = vpot; + + //quadrupolar loops + Acov[3]=vpot*sin((M_PI/2.-th_mag)); + +#elif (NTORUS==14) + //torus from a=0 harm init.c + dipolar oops + ldouble lambda = 15.; + ldouble anorm=1.; //BOBMARK: not used, letting HARM normalize the field + ldouble rchop = 550.; //outer boundary of field loops + ldouble u_av = pp[UU]; + ldouble u_av_chop, u_av_mid; + //midplane at r_mag + init_dsandvels_limotorus(r_mag, M_PI/2., BHSPIN, &rho, &u_av_mid, &ell); + //midplane at rchop + init_dsandvels_limotorus(rchop, M_PI/2., BHSPIN, &rho, &u_av_chop, &ell); + + //vector potential follows contours of UU + ldouble uchop = u_av - u_av_chop; //vpot->zero on contour of radius r=rchop + ldouble uchopmid = u_av_mid - u_av_chop; //vpot->zero away from midplane + + ldouble rin=LT_RIN; + ldouble STARTFIELD = 2.*rin; + ldouble q, fr, fr_start, vpot=0.; + if (r_mag > STARTFIELD && r_mag < rchop) { + q = anorm * (uchop - 0.2*uchopmid) / (0.8*uchopmid) * pow(sin(th_mag), 3); + } else q = 0; + + if(q > 0.) { + fr = (pow(r_mag,0.6)/0.6 + 0.5/0.4*pow(r_mag,-0.4)) / lambda; + fr_start = (pow(STARTFIELD,0.6)/0.6 + 0.5/0.4*pow(STARTFIELD,-0.4)) / lambda; + vpot += q * sin(fr - fr_start) ; + } + + Acov[3]=vpot; + +#else //standard single poloidal loop + ldouble dmax = rhoCGS2GU(RHOMAX_VPOT_CGS)*RMAX_VPOT*RMAX_VPOT; + init_dsandvels_limotorus(r_mag, th_mag, BHSPIN, &rho, &uint, &ell); + if(fabs((rho-pp[RHO])/pp[RHO])>.5) rho=pp[RHO]; + Acov[3]=my_max(pow(rho*r_mag*r_mag/dmax , 2.) - RHO_VPOT_CUT, 0.) * pow(sin(fabs(th_mag)), 4.); +#endif + pp[B1]=Acov[1]; + pp[B2]=Acov[2]; + pp[B3]=Acov[3]; +#endif + +} //if rho<0 + +//entropy +pp[ENTR]=calc_Sfromu(pp[RHO],pp[UU],ix,iy,iz); + +//thermodynamics and relativistic electrons +#ifdef EVOLVEELECTRONS +#ifdef RELELECTRONS +int ie; +for (ie=0; ie < NRELBIN; ie++) +{ + if(relel_gammas[ie]RELEL_INJ_MAX) + pp[NEREL(ie)] = 0.0 ; //Always initialize with zero nonthermal + else + pp[NEREL(ie)] = pow(relel_gammas[ie],-RELEL_HEAT_INDEX); +} + +ldouble unth_tmp = calc_relel_uint(pp); +if (unth_tmp > 0) +{ + ldouble scalefac = 1.e-5*pp[UU]/unth_tmp; + for(ie=0;iemaxbeta) + { + maxbeta=pgas/pmag; + //printf("%d %d > %e %e %e\n",ix,iy,pmag,pgas,maxbeta); + } + if(pmag>maxpmag) + { + maxpmag=pmag; + //printf("%d %d > %e %e %e\n",ix,iy,pmag,pgas,maxbeta); + } + if(pgas>maxpgas) + { + maxpgas=pgas; + //printf("%d %d > %e %e %e\n",ix,iy,pmag,pgas,maxbeta); + } + +#else //normalizing wrt to the equatorial plane + if(geom.iy==NY/2) + { +#pragma omp critical + if(pgas/pmag>maxbeta) + { + maxbeta=pgas/pmag; + //printf("%d %d > %e %e %e\n",ix,iy,pmag,pgas,maxbeta); + } + if(pmag>maxpmag) + { + maxpmag=pmag; + //printf("%d %d > %e %e %e\n",ix,iy,pmag,pgas,maxbeta); + } + if(pgas>maxpgas) + { + maxpgas=pgas; + //printf("%d %d > %e %e %e\n",ix,iy,pmag,pgas,maxbeta); + } + + } +#endif + } + } + } + + ldouble global_maxbeta=maxbeta; + ldouble global_maxpgas=maxpgas; + ldouble global_maxpmag=maxpmag; + +//choosing maximal maxbeta from the whole domain +#ifdef MPI + MPI_Barrier(MPI_COMM_WORLD); + + MPI_Allreduce(&maxbeta, &global_maxbeta, 1, MPI_DOUBLE, MPI_MAX, + MPI_COMM_WORLD); + MPI_Allreduce(&maxpgas, &global_maxpgas, 1, MPI_DOUBLE, MPI_MAX, + MPI_COMM_WORLD); + MPI_Allreduce(&maxpmag, &global_maxpmag, 1, MPI_DOUBLE, MPI_MAX, + MPI_COMM_WORLD); + + //maxbeta=global_maxbeta; +#endif + +maxbeta = global_maxpgas/global_maxpmag; +ldouble fac=sqrt(maxbeta/MAXBETA); + +//manual normalization - verify beta! +#ifdef BETANORMFACTOR +fac=BETANORMFACTOR; +#endif +//printf("rescaling magn.fields by %e (%e)\n",fac,maxbeta); + +#pragma omp parallel for private(ix,iy,iz) schedule (dynamic) +for(iz=0;iz=rbreak2 ? rbreak2 : r , a) ); +} + +FTYPE rtbis(FTYPE (*func)(FTYPE,FTYPE*), FTYPE *parms, FTYPE x1, FTYPE x2, FTYPE xacc) +//Taken from HARM:nrutil.c +{ + int j; + FTYPE dx,f,fmid,xmid,rtb; + f=(*func)(x1, parms); + fmid=(*func)(x2, parms); + if (f*fmid >= 0.0) { + printf( "f(%g)=%g f(%g)=%g\n", x1, f, x2, fmid ); + printf("Root must be bracketed for bisection in rtbis\n"); + } + rtb = (f < 0.0) ? (dx=x2-x1,x1) : (dx=x1-x2,x2); //Orient the search so that f>0 lies at x+dx. + for (j=1;j<=100;j++) { + fmid=(*func)(xmid=rtb+(dx *= 0.5),parms); //Bisection loop. + if (fmid <= 0.0) { + rtb=xmid; + } + if (fabs(dx) < xacc || fmid == 0.0) { + return rtb; + } + } + printf("Too many bisections in rtbis\n"); + return 0.0; //Never get here. +} + + +// von Zeipel cylinder radius corresponding to equatorial radius +FTYPE Lam_of_rEq(FTYPE req, FTYPE a, FTYPE rbreak1, FTYPE rbreak2, FTYPE xi) { + + FTYPE gdtt,gdtp,gdpp,l,lam,lam2; + compute_gd(req, M_PI_2, a, &gdtt, &gdtp, &gdpp); + l = l3d(req, a, rbreak1, rbreak2, xi); + lam2 = -l*(l*gdtp+gdpp)/(l*gdtt+gdtp); + lam = sqrt(lam2); + + return lam; + +} +FTYPE rEqofLam_func(FTYPE req, FTYPE *parms) { + FTYPE lamguess, lam, a, rbreak1, rbreak2, xi; + + lam = parms[0]; + a = parms[1]; + rbreak1 = parms[2]; + rbreak2 = parms[3]; + xi = parms[4]; + + lamguess = Lam_of_rEq(req, a, rbreak1, rbreak2, xi); + + return ( lam*lam - lamguess*lamguess); +} + +// equatorial radius corresponding to von Zeipel cylinder radius +FTYPE rEq_of_Lam(FTYPE lam, FTYPE a, FTYPE rbreak1, FTYPE rbreak2, FTYPE xi) { + + FTYPE parms[5]; + + //store params in an array before function call + parms[0] = lam; + parms[1] = a; + parms[2] = rbreak1; + parms[3] = rbreak2; + parms[4] = xi; + + //solve for equatorial radius using bisection, specify large enough root search range, (1e-3, 1e3) + //demand accuracy 5x machine prec. + //in non-rel limit req = lambda , use 10x that as the upper limit: + + + if(rEqofLam_func(.01*lam,parms)*rEqofLam_func(10*lam,parms)>0) + return -1.; + // printf("%e %e %e %e %e \n",lam,Lam_of_rEq(.01*lam, a, rbreak1, rbreak2, xi),Lam_of_rEq(10*lam, a, rbreak1, rbreak2, xi), rEqofLam_func(.01*lam,parms),rEqofLam_func(10*lam,parms)); + else return( rtbis( &rEqofLam_func, parms, .01*lam, 10*lam, 5.*DBL_EPSILON ) ); +} + + + +FTYPE lamBL_func(FTYPE lam, FTYPE *parms) { + FTYPE r,req,th, gdtt, gdtp, gdpp, gdtpeq, a, rbreak1, rbreak2, xi, l; + + gdtt = parms[0]; + gdtp = parms[1]; + gdpp = parms[2]; + a = parms[3]; + rbreak1 = parms[4]; + rbreak2 = parms[5]; + xi = parms[6]; + r = parms[7]; + th = parms[8]; + + req = rEq_of_Lam(lam, a, rbreak1, rbreak2, xi); + l = l3d(req, a, rbreak1, rbreak2, xi); + + return ( lam*lam + l * (l*gdtp + gdpp) / (l*gdtt + gdtp) ); +} + +// von Zeipel cylinder radius of an arbitrary r, theta +FTYPE lamBL(FTYPE r, FTYPE th, FTYPE gdtt, FTYPE gdtp, FTYPE gdpp, FTYPE a, FTYPE rbreak1, FTYPE rbreak2, FTYPE xi) { + FTYPE R = r*sin(th);//, used as initial guess for lamBL + + FTYPE parms[9]; + + //store params in an array before function call + parms[0] = gdtt; + parms[1] = gdtp; + parms[2] = gdpp; + parms[3] = a; + parms[4] = rbreak1; + parms[5] = rbreak2; + parms[6] = xi; + parms[7] = r; + parms[8] = th; + + //solve for lin using bisection, specify large enough root search range, (1e-3, 1e3) + //demand accuracy 5x machine prec. + //in non-rel limit rml = r , use 10x that as the upper limit: + if(th==M_PI_2) + return( Lam_of_rEq(r, a, rbreak1, rbreak2, xi) ); + else + return( rtbis( &lamBL_func, parms, R, 10*R, 5.*DBL_EPSILON ) ); +} + +FTYPE omega3d( FTYPE l, FTYPE gdtt, FTYPE gdtp, FTYPE gdpp ) { + return( -(gdtt*l + gdtp)*pow(gdpp + gdtp*l,-1) ); +} + +FTYPE compute_Agrav( FTYPE om, FTYPE gdtt, FTYPE gdtp, FTYPE gdpp ){ + return (sqrt(fabs(1./ ( gdtt + 2*om*gdtp + pow(om,2)*gdpp ) ))); +} + +/* +FTYPE rmidlam( FTYPE x, FTYPE *parms ) { + FTYPE lamsq, gdtt, gdtp, gdpp, a, rbreak1, rbreak2, xi; + FTYPE lam_x, ans; + + lamsq = parms[0]; // square of target lambda + gdtt = parms[1]; + gdtp = parms[2]; + gdpp = parms[3]; + a = parms[4]; + rbreak1 = parms[5]; + rbreak2 = parms[6]; + xi = parms[7]; + + compute_gd(x, M_PI_2, a, &gdtt, &gdtp, &gdpp); + lam_x = lamBL(x, M_PI_2, gdtt, gdtp, gdpp, a, rbreak1, rbreak2, xi); // lambda at current value of x + + ans = lamsq - lam_x*lam_x; + + return(ans); +} + +FTYPE limotorus_findrml(FTYPE lam, FTYPE gdtt, FTYPE gdtp, FTYPE gdpp, FTYPE a, FTYPE rbreak1, FTYPE rbreak2, FTYPE xi) { + FTYPE parms[8]; + FTYPE rml; + + //store params in an array before function call + parms[0] = lam*lam; + parms[1] = gdtt; + parms[2] = gdtp; + parms[3] = gdpp; + parms[4] = a; + parms[5] = rbreak1; + parms[6] = rbreak2; + parms[7] = xi; + + //solve for lin using bisection, specify large enough root search range, (1e-3, 1e3) + //demand accuracy 5x machine prec. + //in non-rel limit rml = lam , use 10x that as the upper limit: + rml = rtbis( &rmidlam, parms, 3, 10*lam, 5.*DBL_EPSILON ); + + return( rml ); +} +*/ + +struct lnf_int_params { + FTYPE a, rbreak1, rbreak2, xi; +}; + +// Integrand in expression for lnf. This function is called +// by the GSL quadrature routine. +// all quantities are equatorial +FTYPE lnf_integrand(FTYPE r, void *params) { + struct lnf_int_params *pars = (struct lnf_int_params *) params; + FTYPE a = pars->a, rbreak1 = pars->rbreak1, rbreak2 = pars->rbreak2, xi = pars->xi; + FTYPE rhalf, r2, r3, r4,a2, lam, lam2, lamroot; + FTYPE gdtt, gdtp, gdpp; + FTYPE l, om, dl_dr, dx_dr, dom_dr, integrand; + FTYPE term1, term2, term3, term4, D, E; + FTYPE oneplusx, dx_dlam, dlK_dlam, om_numerator, om_denominator; + + // all values below are midplane values + rhalf = sqrt(r); + r2 = r*r; + r3 = r2*r; + r4 = r3*r; + + a2 = a*a; + //compute_gd(r, M_PI_2, a, &gdtt, &gdtp, &gdpp); + //lam = Lam_of_rEq(r, a, rbreak1, rbreak2, xi); + //lam = lamBL(r, M_PI_2, gdtt, gdtp, gdpp, a, rbreak1, rbreak2, xi); + //lam2 = lam*lam; + + l = l3d(r, a, rbreak1, rbreak2, xi); + + term1 = (r-2) * l; + term2 = 2*a; + term3 = r3 + a2*(r+2); + term4 = term2 * l; + + om_numerator = term1 + term2; + om_denominator = term3 - term4; + om = om_numerator / om_denominator; + + //om = omega3d( l, gdtt, gdtp, gdpp ) + + // derivatives + if (r <= rbreak1 || r >= rbreak2) { + dl_dr = 0; + } else { + //dl_dr_A = (a2*r2*rhalf + 4*a*r3 - 3*r3*rhalf +0.5*r4*rhalf); + //dl_dr_B = (a*r - 2*r*rhalf + r2*rhalf); + //dl_dr = xi * dl_dr_A / (dl_dr_B * dl_dr_B); + + oneplusx = 1 - 2/r + a/r/rhalf; + dx_dr = 2/r2 - 1.5*a/r2/rhalf; + dl_dr = (oneplusx*0.5/rhalf + (1-rhalf)*dx_dr) / (oneplusx*oneplusx); + + //old + //oneplusx = 1 - 2/lam + a*pow(lam,-1.5); + //dx_dlam = 2*pow(lam,-2) - 1.5*a*pow(lam,-2.5); + //lamroot = sqrt(lam); + //dlK_dlam = (oneplusx*0.5*pow(lamroot,-1) + (a-lamroot)*dx_dlam) / pow(oneplusx,2); + + //D = term3 - 2*term4 - lam2 * (r-2); + //E = l * (3*r2 + a2 - lam2); + //dl_dr = E / ( 2*lam*om_numerator / (xi*dlK_dlam) - D); + } + + dom_dr = ( om_denominator * (l + (r-2)*dl_dr) - om_numerator * (3*r2+a2+2*a*dl_dr) ) + / pow(om_denominator, 2); + + integrand = -l/(1 - om*l) * dom_dr; + + return( integrand ); +} + +int init_dsandvels_limotorus(FTYPE r, FTYPE th, FTYPE a, FTYPE *rhoout, FTYPE *uuout, FTYPE *ell) +{ + + FTYPE rho=-1.; + FTYPE uu=-1.; + FTYPE l=-1.; + +#ifdef LIMOTORUS + FTYPE kappa; + FTYPE hh, eps; + FTYPE u, ur, uh, up; + int pl; + FTYPE R, rbreak1, rbreak2, xi; + //FTYPE lambreak1, lambreak2; + FTYPE lam, rml, om, Agrav, f3d; + FTYPE lamin, lin, omin, Agravin, f3din; + FTYPE lnf3d, lnferr; + FTYPE gdtt, gdtp, gdpp; + + // GSL quadrature stuff + gsl_integration_workspace *intwork; + int worksize = 1000; // workspace size + gsl_function integrand; + int intstatus; + R = r*sin(th); + if (R < LT_RIN) + { + *rhoout = -1.; + return(0); + } + + /// + /// Parameters + /// + + kappa = LT_KAPPA; // AKMARK: entropy constant that appears in EOS + xi = LT_XI; // AKMARK: omega is set to this fraction of Keplerian omega + rbreak1 = LT_R1; // AKMARK: locations of breaks in torus angular momentum profile + rbreak2 = LT_R2; + + /// + /// Computations at break radii + /// + + // AKMARK: lamBL can calculate lambda at an arbitrary location, but it needs lambreak1,2; + // how to calculate lambreak1,2 themselves? + // Solution: note that lambreak1,2 are only needed by "l3d" in order to determine which region of the torus we are in. + // But lambreak1,2 can be considered to be in region 2, so just need a way to make "l3d" believe that we are in region 2. + // This can be done by setting lambreak1,2 to very small and very large values respectively. + + //compute_gd(rbreak1, M_PI_2, a, &gdtt, &gdtp, &gdpp); + //lambreak1 = lamBL(rbreak1, gdtt, gdtp, gdpp, a, 0, 200000, xi); + //compute_gd(rbreak2, M_PI_2, a, &gdtt, &gdtp, &gdpp); + //lambreak2 = lamBL(rbreak2, gdtt, gdtp, gdpp, a, lambreak1, 200000, xi); + + /// + /// Computations at torus inner edge + /// + + compute_gd(LT_RIN, M_PI_2, a, &gdtt, &gdtp, &gdpp); + lamin = lamBL(LT_RIN, M_PI_2, gdtt, gdtp, gdpp, a, rbreak1, rbreak2, xi); + lin = l3d(LT_RIN, a, rbreak1, rbreak2, xi); + omin = omega3d(lin, gdtt, gdtp, gdpp); + Agravin = compute_Agrav(omin, gdtt, gdtp, gdpp); + f3din = 1.; // the way f3d is defined, it equals 1 at r=rin + + /// + /// Computations at current point: r, th + /// + + compute_gd(r, th, a, &gdtt, &gdtp, &gdpp); + lam = lamBL(r, th, gdtt, gdtp, gdpp, a, rbreak1, rbreak2, xi); + rml = rEq_of_Lam(lam, a, rbreak1, rbreak2, xi); + l = l3d(rml, a, rbreak1, rbreak2, xi); + om = omega3d(l, gdtt, gdtp, gdpp); + Agrav = compute_Agrav(om, gdtt, gdtp, gdpp); + + //rml = limotorus_findrml(lam, gdtt, gdtp, gdpp, a, rbreak1, rbreak2, xi ); + // f3d requires numerical evaluation of an integral + + // First, set up things that GSL quadrature routine requires + struct lnf_int_params pars = {a, rbreak1, rbreak2, xi}; + integrand.function = &lnf_integrand; // integrand: function + integrand.params = &pars; // other parameters to pass to integrand (besides integration variable) + intwork = gsl_integration_workspace_alloc(worksize); // integration workspace + + // then perform integration to obtain ln(f3d) ... + intstatus = gsl_integration_qags(&integrand, LT_RIN, rml, 0., 1.e-8, worksize, intwork, &lnf3d, &lnferr); + gsl_integration_workspace_free( intwork ); + if (intstatus != GSL_SUCCESS) { + printf("GSL integration failed during limotorus setup at r=%21.15g, th=%21.15g; setting density to zero at this point", r, th); + lnf3d = GSL_NEGINF; // cause density to be 0 at the current point + } + + // ... and finally, f3d + f3d = exp(-lnf3d); + + hh = f3din*Agrav / (f3d*Agravin); + eps = (-1 + hh)*pow(LT_GAMMA,-1); + + if (eps < 0) + rho = -1; + else + rho = pow((-1 + LT_GAMMA)*eps*pow(kappa,-1),pow(-1 + LT_GAMMA,-1)); + + uu = kappa * pow(rho, LT_GAMMA) / (LT_GAMMA - 1.); +#endif + + *rhoout = rho; + *ell=l; + *uuout = uu; + + //my_err("wfsg"); + //printf("%e %e %e\n",rho,kappa,*uuout); + //exit(0); getch(); + + return(0); +} + +// Calculate approximate density maximum for limotorus +// Using penna et al eq 24 +FTYPE rhomax_limotorus_approx(FTYPE a) +{ + ldouble rho=-1.; +#ifdef LIMOTORUS + ldouble ell1=lK(LT_R1, a); // angular momentum at transition point + ldouble rmax1 = LT_XI*LT_XI*ell1*ell1 - 4.; // location of pressure/density maximum + ldouble rmax; + if(rmax1 < LT_RIN) + rmax = LT_RIN; + else + rmax = rmax1; + + ldouble uint,ell; + init_dsandvels_limotorus(rmax, 0.5*M_PI, a, &rho, &uint, &ell); // density at the maximum + //printf("rmax %.4f %.4f %.2e \n", rmax1, rmax, rho); +#endif + return rho; +} + + +// Fishbone & Moncrief (ApJ, 207, 962, 1976) torus with inner edge at rin and pressure maximum at rmax. + +int init_dsandvels_fishbone_moncrief(FTYPE r, FTYPE th, FTYPE a, FTYPE *rhoout, FTYPE *uuout, FTYPE *ellout) +{ + // FM_rin, FM_rmax are the equatorial radii of the torus inner endge and pressure maximum + ldouble rho=-1.; + ldouble uu=-1.; + ldouble ell=-1.; + +#ifndef LIMOTORUS + ldouble rin = FM_rin; + ldouble rmax = FM_rmax; + + ldouble Delta = r*r - 2.*r + a*a; + ldouble Sigma = r*r + a*a*cos(th)*cos(th); + ldouble A = (r*r + a*a)*(r*r + a*a) - Delta*a*a*sin(th)*sin(th); + + // The constant l is evalulated at rmax using eq (3.8) in FM76 + + ell = (rmax*rmax*rmax*rmax + rmax*rmax*a*a - 2.*rmax*a*a - a*sqrt(rmax)*(rmax*rmax - a*a)) / (rmax*rmax - 3.*rmax + 2.*a*sqrt(rmax)) / pow(rmax,1.5); + + //ln h(r,th) (eq (3.6) in FM76) + + ldouble hlog = 0.5*log((1. + sqrt(1. + 4.*ell*ell*Sigma*Sigma*Delta/(A*sin(th)*A*sin(th))))/(Sigma*Delta/A)) - 0.5*sqrt(1 + 4.*ell*ell*Sigma*Sigma*Delta/(A*sin(th)*A*sin(th))) - 2.*a*r*ell/A - 0.5*log((1. + sqrt(1. + 4.*ell*ell*rin*rin*(rin*rin - 2.*rin + a*a)*pow((rin*rin*rin + rin*a*a + 2.*a*a),-2)))/(rin*(rin*rin - 2.*rin + a*a)/(rin*rin*rin + rin*a*a + 2.*a*a))) + 0.5*sqrt(1. + 4.*ell*ell*rin*rin*(rin*rin - 2.*rin + a*a)/pow((rin*rin*rin + rin*a*a + 2.*a*a),2)) + ell*2.*a/(rin*rin*rin + rin*a*a + 2.*a*a); + ldouble h = exp(hlog); + + ldouble Delta0 = rmax*rmax - 2.*rmax + a*a; + ldouble Sigma0 = rmax*rmax; + ldouble A0 = (rmax*rmax + a*a)*(rmax*rmax + a*a) - Delta0*a*a; + + //ln h(rmax,pi/2) (eq (3.6) in FM76) + + ldouble hlog0 = 0.5*log((1. + sqrt(1. + 4.*ell*ell*Sigma0*Sigma0*Delta0/(A0*A0)))/(Sigma0*Delta0/A0)) - 0.5*sqrt(1 + 4.*ell*ell*Sigma0*Sigma0*Delta0/(A0*A0)) - 2.*a*rmax*ell/A0 - 0.5*log((1. + sqrt(1. + 4.*ell*ell*rin*rin*(rin*rin - 2.*rin + a*a)*pow((rin*rin*rin + rin*a*a + 2.*a*a),-2)))/(rin*(rin*rin - 2.*rin + a*a)/(rin*rin*rin + rin*a*a + 2.*a*a))) + 0.5*sqrt(1. + 4.*ell*ell*rin*rin*(rin*rin - 2.*rin + a*a)/pow((rin*rin*rin + rin*a*a + 2.*a*a),2)) + ell*2.*a/(rin*rin*rin + rin*a*a + 2.*a*a); + ldouble h0 = exp(hlog0); + + ldouble n = 1. / (GAMMA - 1.); + ldouble p0, p; + if (h >= 1. && r >= rin) // the point is inside the torus + { + // Set rho = FM_rho0 ((h-1)/(h0-1))^n (note: FM_rho0 is the maximum density) + // p0 = FM_rho0 (h0-1)/(n+1) (note: p0 is the maximum pressure) + // p = p0 ((h-1)/(h0-1))^(n+1) + // u = p/(Gamma-1) + + rho = FM_rho0 * pow(((h - 1.) / (h0 - 1.)),n); + p0 = FM_rho0 * (h0 - 1.) / (n + 1.); + p = p0 * pow(((h - 1.) / (h0 - 1.)),(n+1)); + uu = p / (GAMMA - 1.); + } + else // the point is outside the torus + { + // Set rho = u = -1 to indicate to the calling program that this is outside the torus + + rho = -1.; + uu = -1.; + } + + // Copy and return results +#endif + *rhoout = rho; + *uuout = uu; + *ellout = ell; + + //printf("%e %e %e %e %e %e\n", r, h-1., h0-1., rho, p, ell); + + return 0; +} diff --git a/PROBLEMS/RADSURVEY/tools_fish.c b/PROBLEMS/RADSURVEY/tools_fish.c new file mode 100755 index 0000000..bdba06d --- /dev/null +++ b/PROBLEMS/RADSURVEY/tools_fish.c @@ -0,0 +1,64 @@ +// Fishbone & Moncrief (ApJ, 207, 962, 1976) torus with inner edge at rin and pressure maximum at rmax. + +int init_dsandvels_fishbone_moncrief(FTYPE r, FTYPE th, FTYPE a, FTYPE *rhoout, FTYPE *uuout, FTYPE *ellout) +{ + // FM_rin, FM_rmax are the equatorial radii of the torus inner endge and pressure maximum + + ldouble rin = FM_rin; + ldouble rmax = FM_rmax; + + ldouble Delta = r*r - 2.*r + a*a; + ldouble Sigma = r*r + a*a*cos(th)*cos(th); + ldouble A = (r*r + a*a)*(r*r + a*a) - Delta*a*a*sin(th)*sin(th); + + // The constant l is evalulated at rmax using eq (3.8) in FM76 + + ldouble ell = (rmax*rmax*rmax*rmax + rmax*rmax*a*a - 2.*rmax*a*a - a*sqrt(rmax)*(rmax*rmax - a*a)) / (rmax*rmax - 3.*rmax + 2.*a*sqrt(rmax)) / pow(rmax,1.5); + + //ln h(r,th) (eq (3.6) in FM76) + + ldouble hlog = 0.5*log((1. + sqrt(1. + 4.*ell*ell*Sigma*Sigma*Delta/(A*sin(th)*A*sin(th))))/(Sigma*Delta/A)) - 0.5*sqrt(1 + 4.*ell*ell*Sigma*Sigma*Delta/(A*sin(th)*A*sin(th))) - 2.*a*r*ell/A - 0.5*log((1. + sqrt(1. + 4.*ell*ell*rin*rin*(rin*rin - 2.*rin + a*a)*pow((rin*rin*rin + rin*a*a + 2.*a*a),-2)))/(rin*(rin*rin - 2.*rin + a*a)/(rin*rin*rin + rin*a*a + 2.*a*a))) + 0.5*sqrt(1. + 4.*ell*ell*rin*rin*(rin*rin - 2.*rin + a*a)/pow((rin*rin*rin + rin*a*a + 2.*a*a),2)) + ell*2.*a/(rin*rin*rin + rin*a*a + 2.*a*a); + ldouble h = exp(hlog); + + ldouble Delta0 = rmax*rmax - 2.*rmax + a*a; + ldouble Sigma0 = rmax*rmax; + ldouble A0 = (rmax*rmax + a*a)*(rmax*rmax + a*a) - Delta0*a*a; + + //ln h(rmax,pi/2) (eq (3.6) in FM76) + + ldouble hlog0 = 0.5*log((1. + sqrt(1. + 4.*ell*ell*Sigma0*Sigma0*Delta0/(A0*A0)))/(Sigma0*Delta0/A0)) - 0.5*sqrt(1 + 4.*ell*ell*Sigma0*Sigma0*Delta0/(A0*A0)) - 2.*a*rmax*ell/A0 - 0.5*log((1. + sqrt(1. + 4.*ell*ell*rin*rin*(rin*rin - 2.*rin + a*a)*pow((rin*rin*rin + rin*a*a + 2.*a*a),-2)))/(rin*(rin*rin - 2.*rin + a*a)/(rin*rin*rin + rin*a*a + 2.*a*a))) + 0.5*sqrt(1. + 4.*ell*ell*rin*rin*(rin*rin - 2.*rin + a*a)/pow((rin*rin*rin + rin*a*a + 2.*a*a),2)) + ell*2.*a/(rin*rin*rin + rin*a*a + 2.*a*a); + ldouble h0 = exp(hlog0); + + ldouble n = 1. / (GAMMA - 1.); + ldouble rho, uu, p0, p; + if (h >= 1. && r >= rin) // the point is inside the torus + { + // Set rho = FM_rho0 ((h-1)/(h0-1))^n (note: FM_rho0 is the maximum density) + // p0 = FM_rho0 (h0-1)/(n+1) (note: p0 is the maximum pressure) + // p = p0 ((h-1)/(h0-1))^(n+1) + // u = p/(Gamma-1) + + rho = FM_rho0 * pow(((h - 1.) / (h0 - 1.)),n); + p0 = FM_rho0 * (h0 - 1.) / (n + 1.); + p = p0 * pow(((h - 1.) / (h0 - 1.)),(n+1)); + uu = p / (GAMMA - 1.); + } + else // the point is outside the torus + { + // Set rho = u = -1 to indicate to the calling program that this is outside the torus + + rho = -1.; + uu = -1.; + } + + // Copy and return results + + *rhoout = rho; + *uuout = uu; + *ellout = ell; + + //printf("%e %e %e %e %e %e\n", r, h-1., h0-1., rho, p, ell); + + return 0; +} + diff --git a/PROBLEMS/RADSURVEY/tools_limo.c b/PROBLEMS/RADSURVEY/tools_limo.c new file mode 100755 index 0000000..6eebfe3 --- /dev/null +++ b/PROBLEMS/RADSURVEY/tools_limo.c @@ -0,0 +1,443 @@ +void compute_gd( FTYPE r, FTYPE th, FTYPE a, FTYPE *gdtt, FTYPE *gdtp, FTYPE *gdpp ) { + FTYPE Sigma, tmp; + + Sigma = (pow(a*cos(th),2) + pow(r,2)); + tmp = 2*r*pow(sin(th),2)/Sigma; + + //metric (expressions taken from limotorus4.nb): + *gdtt = -1 + 2*r/Sigma; + *gdtp = -a*tmp; + *gdpp = (pow(r,2) + pow(a,2)*(1+tmp)) * pow(sin(th),2); +} + +// Keplerian, equatorial angular momentum density: l = u_phi / u_t +FTYPE lK(FTYPE r, FTYPE a) { + FTYPE curlyF, curlyG; + + curlyF = 1 - 2*a/pow(r, 1.5) + pow(a/r, 2); + curlyG = 1 - 2/r + a/pow(r, 1.5); + return( pow(r,0.5) * curlyF / curlyG ); +} + +FTYPE l3d(FTYPE r, FTYPE a, FTYPE rbreak1, FTYPE rbreak2, FTYPE xi) { + return ( xi * lK( r<=rbreak1 ? rbreak1 : r>=rbreak2 ? rbreak2 : r , a) ); +} + +//FTYPE l3d(FTYPE lam, FTYPE a, FTYPE lambreak1, FTYPE lambreak2, FTYPE xi) { +// return ( xi * lK( lam<=lambreak1 ? lambreak1 : lam>=lambreak2 ? lambreak2 : lam , a) ); +//} + +FTYPE rtbis(FTYPE (*func)(FTYPE,FTYPE*), FTYPE *parms, FTYPE x1, FTYPE x2, FTYPE xacc) +//Taken from HARM:nrutil.c +{ + int j; + FTYPE dx,f,fmid,xmid,rtb; + f=(*func)(x1, parms); + fmid=(*func)(x2, parms); + if (f*fmid >= 0.0) { + printf( "f(%g)=%g f(%g)=%g\n", x1, f, x2, fmid ); + printf("Root must be bracketed for bisection in rtbis\n"); + } + rtb = (f < 0.0) ? (dx=x2-x1,x1) : (dx=x1-x2,x2); //Orient the search so that f>0 lies at x+dx. + for (j=1;j<=100;j++) { + fmid=(*func)(xmid=rtb+(dx *= 0.5),parms); //Bisection loop. + if (fmid <= 0.0) { + rtb=xmid; + } + if (fabs(dx) < xacc || fmid == 0.0) { + return rtb; + } + } + printf("Too many bisections in rtbis\n"); + return 0.0; //Never get here. +} + + +// von Zeipel cylinder radius corresponding to equatorial radius +FTYPE Lam_of_rEq(FTYPE req, FTYPE a, FTYPE rbreak1, FTYPE rbreak2, FTYPE xi) { + + FTYPE gdtt,gdtp,gdpp,l,lam,lam2; + compute_gd(req, M_PI_2, a, &gdtt, &gdtp, &gdpp); + l = l3d(req, a, rbreak1, rbreak2, xi); + lam2 = -l*(l*gdtp+gdpp)/(l*gdtt+gdtp); + lam = sqrt(lam2); + + return lam; + +} +FTYPE rEqofLam_func(FTYPE req, FTYPE *parms) { + FTYPE lamguess, lam, a, rbreak1, rbreak2, xi; + + lam = parms[0]; + a = parms[1]; + rbreak1 = parms[2]; + rbreak2 = parms[3]; + xi = parms[4]; + + lamguess = Lam_of_rEq(req, a, rbreak1, rbreak2, xi); + + return ( lam*lam - lamguess*lamguess); +} + +// equatorial radius corresponding to von Zeipel cylinder radius +FTYPE rEq_of_Lam(FTYPE lam, FTYPE a, FTYPE rbreak1, FTYPE rbreak2, FTYPE xi) { + + FTYPE parms[5]; + + //store params in an array before function call + parms[0] = lam; + parms[1] = a; + parms[2] = rbreak1; + parms[3] = rbreak2; + parms[4] = xi; + + //solve for equatorial radius using bisection, specify large enough root search range, (1e-3, 1e3) + //demand accuracy 5x machine prec. + //in non-rel limit req = lambda , use 10x that as the upper limit: + + + if(rEqofLam_func(.01*lam,parms)*rEqofLam_func(10*lam,parms)>0) + return -1.; + // printf("%e %e %e %e %e \n",lam,Lam_of_rEq(.01*lam, a, rbreak1, rbreak2, xi),Lam_of_rEq(10*lam, a, rbreak1, rbreak2, xi), rEqofLam_func(.01*lam,parms),rEqofLam_func(10*lam,parms)); + else return( rtbis( &rEqofLam_func, parms, .01*lam, 10*lam, 5.*DBL_EPSILON ) ); +} + + + +FTYPE lamBL_func(FTYPE lam, FTYPE *parms) { + FTYPE r,req,th, gdtt, gdtp, gdpp, gdtpeq, a, rbreak1, rbreak2, xi, l; + + gdtt = parms[0]; + gdtp = parms[1]; + gdpp = parms[2]; + a = parms[3]; + rbreak1 = parms[4]; + rbreak2 = parms[5]; + xi = parms[6]; + r = parms[7]; + th = parms[8]; + + req = rEq_of_Lam(lam, a, rbreak1, rbreak2, xi); + l = l3d(req, a, rbreak1, rbreak2, xi); + + return ( lam*lam + l * (l*gdtp + gdpp) / (l*gdtt + gdtp) ); +} + +// von Zeipel cylinder radius of an arbitrary r, theta +FTYPE lamBL(FTYPE r, FTYPE th, FTYPE gdtt, FTYPE gdtp, FTYPE gdpp, FTYPE a, FTYPE rbreak1, FTYPE rbreak2, FTYPE xi) { + FTYPE R = r*sin(th);//, used as initial guess for lamBL + + FTYPE parms[9]; + + //store params in an array before function call + parms[0] = gdtt; + parms[1] = gdtp; + parms[2] = gdpp; + parms[3] = a; + parms[4] = rbreak1; + parms[5] = rbreak2; + parms[6] = xi; + parms[7] = r; + parms[8] = th; + + //solve for lin using bisection, specify large enough root search range, (1e-3, 1e3) + //demand accuracy 5x machine prec. + //in non-rel limit rml = r , use 10x that as the upper limit: + if(th==M_PI_2) + return( Lam_of_rEq(r, a, rbreak1, rbreak2, xi) ); + else + return( rtbis( &lamBL_func, parms, R, 10*R, 5.*DBL_EPSILON ) ); +} + +FTYPE omega3d( FTYPE l, FTYPE gdtt, FTYPE gdtp, FTYPE gdpp ) { + return( -(gdtt*l + gdtp)*pow(gdpp + gdtp*l,-1) ); +} + +FTYPE compute_Agrav( FTYPE om, FTYPE gdtt, FTYPE gdtp, FTYPE gdpp ){ + return (sqrt(fabs(1./ ( gdtt + 2*om*gdtp + pow(om,2)*gdpp ) ))); +} + +/* +FTYPE rmidlam( FTYPE x, FTYPE *parms ) { + FTYPE lamsq, gdtt, gdtp, gdpp, a, rbreak1, rbreak2, xi; + FTYPE lam_x, ans; + + lamsq = parms[0]; // square of target lambda + gdtt = parms[1]; + gdtp = parms[2]; + gdpp = parms[3]; + a = parms[4]; + rbreak1 = parms[5]; + rbreak2 = parms[6]; + xi = parms[7]; + + compute_gd(x, M_PI_2, a, &gdtt, &gdtp, &gdpp); + lam_x = lamBL(x, M_PI_2, gdtt, gdtp, gdpp, a, rbreak1, rbreak2, xi); // lambda at current value of x + + ans = lamsq - lam_x*lam_x; + + return(ans); +} + +FTYPE limotorus_findrml(FTYPE lam, FTYPE gdtt, FTYPE gdtp, FTYPE gdpp, FTYPE a, FTYPE rbreak1, FTYPE rbreak2, FTYPE xi) { + FTYPE parms[8]; + FTYPE rml; + + //store params in an array before function call + parms[0] = lam*lam; + parms[1] = gdtt; + parms[2] = gdtp; + parms[3] = gdpp; + parms[4] = a; + parms[5] = rbreak1; + parms[6] = rbreak2; + parms[7] = xi; + + //solve for lin using bisection, specify large enough root search range, (1e-3, 1e3) + //demand accuracy 5x machine prec. + //in non-rel limit rml = lam , use 10x that as the upper limit: + rml = rtbis( &rmidlam, parms, 3, 10*lam, 5.*DBL_EPSILON ); + + return( rml ); +} +*/ + +struct lnf_int_params { + FTYPE a, rbreak1, rbreak2, xi; +}; + +// Integrand in expression for lnf. This function is called +// by the GSL quadrature routine. +// all quantities are equatorial +FTYPE lnf_integrand(FTYPE r, void *params) { + struct lnf_int_params *pars = (struct lnf_int_params *) params; + FTYPE a = pars->a, rbreak1 = pars->rbreak1, rbreak2 = pars->rbreak2, xi = pars->xi; + FTYPE rhalf, r2, r3, r4,a2, lam, lam2, lamroot; + FTYPE gdtt, gdtp, gdpp; + FTYPE l, om, dl_dr, dx_dr, dom_dr, integrand; + FTYPE term1, term2, term3, term4, D, E; + FTYPE oneplusx, dx_dlam, dlK_dlam, om_numerator, om_denominator; + + // all values below are midplane values + rhalf = sqrt(r); + r2 = r*r; + r3 = r2*r; + r4 = r3*r; + + a2 = a*a; + //compute_gd(r, M_PI_2, a, &gdtt, &gdtp, &gdpp); + //lam = Lam_of_rEq(r, a, rbreak1, rbreak2, xi); + //lam = lamBL(r, M_PI_2, gdtt, gdtp, gdpp, a, rbreak1, rbreak2, xi); + //lam2 = lam*lam; + + l = l3d(r, a, rbreak1, rbreak2, xi); + + term1 = (r-2) * l; + term2 = 2*a; + term3 = r3 + a2*(r+2); + term4 = term2 * l; + + om_numerator = term1 + term2; + om_denominator = term3 - term4; + om = om_numerator / om_denominator; + + //om = omega3d( l, gdtt, gdtp, gdpp ) + + // derivatives + if (r <= rbreak1 || r >= rbreak2) { + dl_dr = 0; + } else { + //dl_dr_A = (a2*r2*rhalf + 4*a*r3 - 3*r3*rhalf +0.5*r4*rhalf); + //dl_dr_B = (a*r - 2*r*rhalf + r2*rhalf); + //dl_dr = xi * dl_dr_A / (dl_dr_B * dl_dr_B); + + oneplusx = 1 - 2/r + a/r/rhalf; + dx_dr = 2/r2 - 1.5*a/r2/rhalf; + dl_dr = (oneplusx*0.5/rhalf + (1-rhalf)*dx_dr) / (oneplusx*oneplusx); + + //old + //oneplusx = 1 - 2/lam + a*pow(lam,-1.5); + //dx_dlam = 2*pow(lam,-2) - 1.5*a*pow(lam,-2.5); + //lamroot = sqrt(lam); + //dlK_dlam = (oneplusx*0.5*pow(lamroot,-1) + (a-lamroot)*dx_dlam) / pow(oneplusx,2); + + //D = term3 - 2*term4 - lam2 * (r-2); + //E = l * (3*r2 + a2 - lam2); + //dl_dr = E / ( 2*lam*om_numerator / (xi*dlK_dlam) - D); + } + + dom_dr = ( om_denominator * (l + (r-2)*dl_dr) - om_numerator * (3*r2+a2+2*a*dl_dr) ) + / pow(om_denominator, 2); + + integrand = -l/(1 - om*l) * dom_dr; + + return( integrand ); +} + +int init_dsandvels_limotorus(FTYPE r, FTYPE th, FTYPE a, FTYPE *rhoout, FTYPE *uuout, FTYPE *ell) +{ + FTYPE kappa; + FTYPE hh, eps; + FTYPE rho, u, ur, uh, up; + int pl; + FTYPE R, rbreak1, rbreak2, xi; + //FTYPE lambreak1, lambreak2; + FTYPE lam, rml, l, om, Agrav, f3d; + FTYPE lamin, lin, omin, Agravin, f3din; + FTYPE lnf3d, lnferr; + FTYPE gdtt, gdtp, gdpp; + + // GSL quadrature stuff + gsl_integration_workspace *intwork; + int worksize = 1000; // workspace size + gsl_function integrand; + int intstatus; + R = r*sin(th); + if (R < LT_RIN) {*rhoout = -1.; return(0);} + + /// + /// Parameters + /// + + kappa = LT_KAPPA; // AKMARK: entropy constant that appears in EOS + xi = LT_XI; // AKMARK: omega is set to this fraction of Keplerian omega + rbreak1 = LT_R1; // AKMARK: locations of breaks in torus angular momentum profile + rbreak2 = LT_R2; + + /// + /// Computations at break radii + /// + + // AKMARK: lamBL can calculate lambda at an arbitrary location, but it needs lambreak1,2; + // how to calculate lambreak1,2 themselves? + // Solution: note that lambreak1,2 are only needed by "l3d" in order to determine which region of the torus we are in. + // But lambreak1,2 can be considered to be in region 2, so just need a way to make "l3d" believe that we are in region 2. + // This can be done by setting lambreak1,2 to very small and very large values respectively. + + //compute_gd(rbreak1, M_PI_2, a, &gdtt, &gdtp, &gdpp); + //lambreak1 = lamBL(rbreak1, gdtt, gdtp, gdpp, a, 0, 200000, xi); + //compute_gd(rbreak2, M_PI_2, a, &gdtt, &gdtp, &gdpp); + //lambreak2 = lamBL(rbreak2, gdtt, gdtp, gdpp, a, lambreak1, 200000, xi); + + /// + /// Computations at torus inner edge + /// + + compute_gd(LT_RIN, M_PI_2, a, &gdtt, &gdtp, &gdpp); + lamin = lamBL(LT_RIN, M_PI_2, gdtt, gdtp, gdpp, a, rbreak1, rbreak2, xi); + lin = l3d(LT_RIN, a, rbreak1, rbreak2, xi); + omin = omega3d(lin, gdtt, gdtp, gdpp); + Agravin = compute_Agrav(omin, gdtt, gdtp, gdpp); + f3din = 1.; // the way f3d is defined, it equals 1 at r=rin + + /// + /// Computations at current point: r, th + /// + + compute_gd(r, th, a, &gdtt, &gdtp, &gdpp); + lam = lamBL(r, th, gdtt, gdtp, gdpp, a, rbreak1, rbreak2, xi); + rml = rEq_of_Lam(lam, a, rbreak1, rbreak2, xi); + l = l3d(rml, a, rbreak1, rbreak2, xi); + om = omega3d(l, gdtt, gdtp, gdpp); + Agrav = compute_Agrav(om, gdtt, gdtp, gdpp); + + //rml = limotorus_findrml(lam, gdtt, gdtp, gdpp, a, rbreak1, rbreak2, xi ); + // f3d requires numerical evaluation of an integral + + // First, set up things that GSL quadrature routine requires + struct lnf_int_params pars = {a, rbreak1, rbreak2, xi}; + integrand.function = &lnf_integrand; // integrand: function + integrand.params = &pars; // other parameters to pass to integrand (besides integration variable) + intwork = gsl_integration_workspace_alloc(worksize); // integration workspace + + // then perform integration to obtain ln(f3d) ... + intstatus = gsl_integration_qags(&integrand, LT_RIN, rml, 0., 1.e-8, worksize, intwork, &lnf3d, &lnferr); + gsl_integration_workspace_free( intwork ); + if (intstatus != GSL_SUCCESS) { + printf("GSL integration failed during limotorus setup at r=%21.15g, th=%21.15g; setting density to zero at this point", r, th); + lnf3d = GSL_NEGINF; // cause density to be 0 at the current point + } + + // ... and finally, f3d + f3d = exp(-lnf3d); + + hh = f3din*Agrav / (f3d*Agravin); + eps = (-1 + hh)*pow(LT_GAMMA,-1); + + if (eps < 0) rho = -1; else rho = pow((-1 + LT_GAMMA)*eps*pow(kappa,-1),pow(-1 + LT_GAMMA,-1)); + + *rhoout = rho; + *ell=l; + *uuout = kappa * pow(rho, LT_GAMMA) / (LT_GAMMA - 1.); + + //my_err("wfsg"); + //printf("%e %e %e\n",rho,kappa,*uuout); + //exit(0); getch(); + + return(0); + +} + +// Fishbone & Moncrief (ApJ, 207, 962, 1976) torus with inner edge at rin and pressure maximum at rmax. + +int init_dsandvels_fishbone_moncrief(FTYPE r, FTYPE th, FTYPE a, FTYPE *rhoout, FTYPE *uuout, FTYPE *ellout) +{ + // FM_rin, FM_rmax are the equatorial radii of the torus inner endge and pressure maximum + + ldouble rin = FM_rin; + ldouble rmax = FM_rmax; + + ldouble Delta = r*r - 2.*r + a*a; + ldouble Sigma = r*r + a*a*cos(th)*cos(th); + ldouble A = (r*r + a*a)*(r*r + a*a) - Delta*a*a*sin(th)*sin(th); + + // The constant l is evalulated at rmax using eq (3.8) in FM76 + + ldouble ell = (rmax*rmax*rmax*rmax + rmax*rmax*a*a - 2.*rmax*a*a - a*sqrt(rmax)*(rmax*rmax - a*a)) / (rmax*rmax - 3.*rmax + 2.*a*sqrt(rmax)) / pow(rmax,1.5); + + //ln h(r,th) (eq (3.6) in FM76) + + ldouble hlog = 0.5*log((1. + sqrt(1. + 4.*ell*ell*Sigma*Sigma*Delta/(A*sin(th)*A*sin(th))))/(Sigma*Delta/A)) - 0.5*sqrt(1 + 4.*ell*ell*Sigma*Sigma*Delta/(A*sin(th)*A*sin(th))) - 2.*a*r*ell/A - 0.5*log((1. + sqrt(1. + 4.*ell*ell*rin*rin*(rin*rin - 2.*rin + a*a)*pow((rin*rin*rin + rin*a*a + 2.*a*a),-2)))/(rin*(rin*rin - 2.*rin + a*a)/(rin*rin*rin + rin*a*a + 2.*a*a))) + 0.5*sqrt(1. + 4.*ell*ell*rin*rin*(rin*rin - 2.*rin + a*a)/pow((rin*rin*rin + rin*a*a + 2.*a*a),2)) + ell*2.*a/(rin*rin*rin + rin*a*a + 2.*a*a); + ldouble h = exp(hlog); + + ldouble Delta0 = rmax*rmax - 2.*rmax + a*a; + ldouble Sigma0 = rmax*rmax; + ldouble A0 = (rmax*rmax + a*a)*(rmax*rmax + a*a) - Delta0*a*a; + + //ln h(rmax,pi/2) (eq (3.6) in FM76) + + ldouble hlog0 = 0.5*log((1. + sqrt(1. + 4.*ell*ell*Sigma0*Sigma0*Delta0/(A0*A0)))/(Sigma0*Delta0/A0)) - 0.5*sqrt(1 + 4.*ell*ell*Sigma0*Sigma0*Delta0/(A0*A0)) - 2.*a*rmax*ell/A0 - 0.5*log((1. + sqrt(1. + 4.*ell*ell*rin*rin*(rin*rin - 2.*rin + a*a)*pow((rin*rin*rin + rin*a*a + 2.*a*a),-2)))/(rin*(rin*rin - 2.*rin + a*a)/(rin*rin*rin + rin*a*a + 2.*a*a))) + 0.5*sqrt(1. + 4.*ell*ell*rin*rin*(rin*rin - 2.*rin + a*a)/pow((rin*rin*rin + rin*a*a + 2.*a*a),2)) + ell*2.*a/(rin*rin*rin + rin*a*a + 2.*a*a); + ldouble h0 = exp(hlog0); + + ldouble n = 1. / (GAMMA - 1.); + ldouble rho, uu, p0, p; + if (h >= 1. && r >= rin) // the point is inside the torus + { + // Set rho = FM_rho0 ((h-1)/(h0-1))^n (note: FM_rho0 is the maximum density) + // p0 = FM_rho0 (h0-1)/(n+1) (note: p0 is the maximum pressure) + // p = p0 ((h-1)/(h0-1))^(n+1) + // u = p/(Gamma-1) + + rho = FM_rho0 * pow(((h - 1.) / (h0 - 1.)),n); + p0 = FM_rho0 * (h0 - 1.) / (n + 1.); + p = p0 * pow(((h - 1.) / (h0 - 1.)),(n+1)); + uu = p / (GAMMA - 1.); + } + else // the point is outside the torus + { + // Set rho = u = -1 to indicate to the calling program that this is outside the torus + + rho = -1.; + uu = -1.; + } + + // Copy and return results + + *rhoout = rho; + *uuout = uu; + *ellout = ell; + + //printf("%e %e %e %e %e %e\n", r, h-1., h0-1., rho, p, ell); + + return 0; +} diff --git a/ko.c b/ko.c index daddc67..c378135 100644 --- a/ko.c +++ b/ko.c @@ -492,7 +492,7 @@ main(int argc, char **argv) free_arrays(); fprint_closefiles(); mpi_myfinalize(); - + return 0; } diff --git a/postproc.c b/postproc.c index 2950f1c..70975e2 100644 --- a/postproc.c +++ b/postproc.c @@ -1301,6 +1301,9 @@ int calc_scalars(ldouble *scalars,ldouble t) ldouble rmri=xxBL[1]/2.; // middle radial cell #if(PROBLEM==69) //INJDISK rmri=20.; +#endif +#if(PROBLEM==RADSURVEY) + rmri=20.; #endif scalars[5]=calc_resmri(rmri); diff --git a/problem.h b/problem.h index dad5dd4..f7ee1d4 100644 --- a/problem.h +++ b/problem.h @@ -140,14 +140,25 @@ //137 KATOTORUS_TILTED - radiative torus initiated like in Kato+2004 with a tilt from BH spin axis //138 JETCOORDS -- test jet coordinate system //139 MADCC -- MAD code comparison +//140 RADSURVEY -- radiative parameter survey //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 107 +#define PROBLEM 140 +#if(PROBLEM==140) + +#define PR_DEFINE "PROBLEMS/RADSURVEY/define.h" +#define PR_BC "PROBLEMS/RADSURVEY/bc.c" +#define PR_INIT "PROBLEMS/RADSURVEY/init.c" +#define PR_KAPPAES "PROBLEMS/RADSURVEY/kappaes.c" +#define PR_TOOLS "PROBLEMS/RADSURVEY/tools.c" +#define PR_POSTINIT "PROBLEMS/RADSURVEY/postinit.c" + +#endif #if(PROBLEM==139) diff --git a/silo.c b/silo.c index 7d3e7f2..63c230a 100644 --- a/silo.c +++ b/silo.c @@ -2,7 +2,6 @@ \brief Routines for SILO output files */ - #ifndef NOSILO //KORAL - silo.c @@ -336,7 +335,7 @@ int fprint_silofile(ldouble time, int num, char* folder, char* prefix) ldouble dxph[3],dx[3]; get_cellsize_out(ix, iy, iz, dx); - /* + /* xx1[0]=0.;xx1[1]=get_xb(ix,0);xx1[2]=get_x(iy,1);xx1[3]=get_x(iz,2); xx2[0]=0.;xx2[1]=get_xb(ix+1,0);xx2[2]=get_x(iy,1);xx2[3]=get_x(iz,2); coco_N(xx1,xx1,MYCOORDS,OUTCOORDS); @@ -352,7 +351,8 @@ int fprint_silofile(ldouble time, int num, char* folder, char* prefix) coco_N(xx1,xx1,MYCOORDS,OUTCOORDS); coco_N(xx2,xx2,MYCOORDS,OUTCOORDS); dx[2]=fabs(xx2[3]-xx1[3]); - */ + */ + dxph[0]=dx[0]*sqrt(geomout.gg[1][1]); dxph[1]=dx[1]*sqrt(geomout.gg[2][2]); dxph[2]=dx[2]*sqrt(geomout.gg[3][3]);