Skip to content

Commit d0aba75

Browse files
pseudotensorpseudotensor
committed
Updated KAZFULL to allow entropy evolution reduction or full entropy evolution. Also introduced new simple table (0.2GB total) that spans very large density and larger temperature range to preserve neutrino physics at lower density (such as neutron decay) that can occur at large radii. Also allows one to avoid reverting to ideal gas EOS with probably generally bad offset and not generally correct adiabatic index.
1 parent 9b6d42e commit d0aba75

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

42 files changed

+948
-445
lines changed

‎coord.c

+3-2
Original file line numberDiff line numberDiff line change
@@ -541,7 +541,7 @@ void write_coord_parms(int defcoordlocal)
541541
else{
542542

543543
// same for all coords (notice no carraige return)
544-
fprintf(out,"%21.15g %21.15g %21.15g %21.15g ",R0,Rin,Rout,hslope);
544+
fprintf(out,"%21.15g %21.15g %21.15g %21.15g %d ",R0,Rin,Rout,hslope,dofull2pi);
545545

546546
if (defcoordlocal == LOGRSINTH) {
547547
}
@@ -636,7 +636,7 @@ void read_coord_parms(int defcoordlocal)
636636
else{
637637
// don't want to overwrite since restart file sets this
638638
// fscanf(in,HEADER4IN,&R0,&Rin,&Rout,&hslope);
639-
fscanf(in,HEADER4IN,&ftemp,&ftemp,&ftemp,&ftemp);
639+
fscanf(in,HEADER5IN,&ftemp,&ftemp,&ftemp,&ftemp,&ftemp);
640640

641641
if (defcoordlocal == LOGRSINTH) {
642642
}
@@ -713,6 +713,7 @@ void read_coord_parms(int defcoordlocal)
713713
MPI_Bcast(&Rin, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
714714
MPI_Bcast(&Rout, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
715715
MPI_Bcast(&hslope, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
716+
MPI_Bcast(&dofull2pi, 1, MPI_INT, MPIid[0], MPI_COMM_GRMHD);
716717

717718
if (defcoordlocal == LOGRSINTH) {
718719
}

‎defs.jon_interp.h

+5
Original file line numberDiff line numberDiff line change
@@ -13,6 +13,11 @@ int numprocs;
1313
int mycpupos[NDIM];
1414
int ncpux1;
1515

16+
// SOME GEOMETRIC VARIABLES (see also global.jon_interp.h for PERIODICINPHI)
17+
int dofull2pi;
18+
int whichdump,whichdumpversion,numcolumns;
19+
20+
1621
int oN1,oN2,oN3,nN1,nN2,nN3 ;
1722
FTYPE refinefactor;
1823
int roN1,roN2,roN3;

‎eos.funcdeclare.h

+7
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,9 @@
11
// some useful wrappers for user simplicity:
22
extern FTYPE pressure_rho0_u_simple(int i, int j, int k, int loc, FTYPE rho, FTYPE u);
3+
extern FTYPE pressure_rho0_u_simple_forcheckinversion(int i, int j, int k, int loc, FTYPE rho, FTYPE u);
34
extern FTYPE cs2_compute_simple(int i, int j, int k, int loc, FTYPE rho, FTYPE u);
45
extern FTYPE compute_entropy_simple(int i, int j, int k, int loc, FTYPE rho, FTYPE u);
6+
extern FTYPE compute_entropy_simple_forcheckinversion(int i, int j, int k, int loc, FTYPE rho, FTYPE u);
57
extern void get_EOS_parms_simple(int*numparms, int i, int j, int k, int loc, FTYPE *parlist);
68
extern void fix_primitive_eos_scalars_simple(int i, int j, int k, int loc, FTYPE *pr);
79

@@ -15,6 +17,7 @@ extern FTYPE dpdu_rho0_u_simple(int i, int j, int k, int loc, FTYPE rho, FTYPE u
1517
// other wrappers
1618
extern int ufromentropy_calc(struct of_geom *ptrgeom, FTYPE entropy, FTYPE *pr);
1719
extern int entropy_calc(struct of_geom *ptrgeom, FTYPE *pr, FTYPE *entropy);
20+
extern int entropy_calc_forcheckinversion(struct of_geom *ptrgeom, FTYPE *pr, FTYPE *entropy);
1821
extern int invertentropyflux_calc(struct of_geom *ptrgeom, FTYPE entropyflux,int dir, struct of_state *q, FTYPE*pr);
1922

2023

@@ -35,6 +38,10 @@ extern FTYPE cs2_compute(int whicheos, FTYPE *EOSextra, FTYPE rho0, FTYPE u);
3538
extern FTYPE compute_entropy(int whicheos, FTYPE *EOSextra, FTYPE rho0, FTYPE u);
3639
extern FTYPE compute_dSdrho(int whicheos, FTYPE *EOSextra, FTYPE rho0, FTYPE u);
3740
extern FTYPE compute_dSdu(int whicheos, FTYPE *EOSextra, FTYPE rho0, FTYPE u);
41+
extern FTYPE compute_specificentropy_wmrho0(int whicheos, FTYPE *EOSextra, FTYPE rho0, FTYPE u);
42+
extern FTYPE compute_dspecificSdrho_wmrho0(int whicheos, FTYPE *EOSextra, FTYPE rho0, FTYPE wmrho0);
43+
extern FTYPE compute_dspecificSdwmrho0_wmrho0(int whicheos, FTYPE *EOSextra, FTYPE rho0, FTYPE wmrho0);
44+
3845
extern FTYPE pressure_wmrho0(int whicheos, FTYPE *EOSextra, FTYPE rho0, FTYPE wmrho0);
3946
extern FTYPE compute_idwmrho0dp(int whicheos, FTYPE *EOSextra, FTYPE rho0, FTYPE wmrho0);
4047
extern FTYPE compute_idrho0dp(int whicheos, FTYPE *EOSextra, FTYPE rho0, FTYPE wmrho0);

‎global.depnmemonics.h

+5-3
Original file line numberDiff line numberDiff line change
@@ -483,13 +483,15 @@ god=deathadflkjasdflkjasdlfkja242424
483483

484484
/* numerical convenience */
485485
#if(REALTYPE==FLOATTYPE)
486-
#define VERYBIG (1.e30)
486+
#define VERYBIG (1.e37)
487+
#define BIG (1.e+35)
488+
#define SMALL (1.e-35)
487489
#else
488490
#define VERYBIG (1.e150)
491+
#define BIG (1.e+300)
492+
#define SMALL (1.e-300)
489493
#endif
490494

491-
#define BIG (1.e+50)
492-
#define SMALL (1.e-50)
493495

494496
#define SLEPSILON (1.e-6)
495497

‎global.jon_interp.h

+5-4
Original file line numberDiff line numberDiff line change
@@ -32,6 +32,7 @@
3232
#define PERIODICINPHI 1 // default: 1
3333

3434

35+
3536
#include <stdio.h>
3637
#include <math.h>
3738
#include <malloc.h>
@@ -166,20 +167,20 @@
166167
// 16 args
167168
// 21 args after going to 3D and doing MBH/QBH
168169
// 6 more args
169-
#define SCANHEADER "%f %d %d %d %f %f %f %f %f %f %f %f %f %f %f %f %f %f %d %f %f %d %d %d %d %d %d"
170+
#define SCANHEADER "%f %d %d %d %f %f %f %f %f %f %f %f %f %f %f %f %f %f %d %f %f %d %d %d %d %d %d %d %d %d"
170171
#elif(REALTYPE==DOUBLETYPE)
171172
#define SCANARG "%lf"
172173
#define SCANARGVEC "%lf %lf %lf"
173174
#define SCANARG4VEC "%lf %lf %lf %lf"
174-
#define SCANHEADER "%lf %d %d %d %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %d %lf %lf %d %d %d %d %d %d"
175+
#define SCANHEADER "%lf %d %d %d %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %d %lf %lf %d %d %d %d %d %d %d %d %d"
175176
#elif(REALTYPE==LONGDOUBLETYPE)
176177
#define SCANARG "%Lf"
177178
#define SCANARGVEC "%Lf %Lf %Lf"
178179
#define SCANARG4VEC "%Lf %Lf %Lf %Lf"
179-
#define SCANHEADER "%Lf %d %d %d %Lf %Lf %Lf %Lf %Lf %Lf %Lf %Lf %Lf %Lf %Lf %Lf %Lf %Lf %d %Lf %Lf %d %d %d %d %d %d "
180+
#define SCANHEADER "%Lf %d %d %d %Lf %Lf %Lf %Lf %Lf %Lf %Lf %Lf %Lf %Lf %Lf %Lf %Lf %Lf %d %Lf %Lf %d %d %d %d %d %d %d %d %d"
180181
#endif
181182

182-
#define PRINTSCANHEADER "%g %d %d %d %g %g %g %g %g %g %g %g %g %g %g %g %g %g %d %g %g %d %d %d %d %d %d\n"
183+
#define PRINTSCANHEADER "%g %d %d %d %g %g %g %g %g %g %g %g %g %g %g %g %g %g %d %g %g %d %d %d %d %d %d %d %d %d\n"
183184

184185

185186

‎init.tools.c

+17-1
Original file line numberDiff line numberDiff line change
@@ -390,7 +390,6 @@ int user1_init_primitives(FTYPE (*prim)[NSTORE2][NSTORE3][NPR], FTYPE (*pstag)[N
390390

391391

392392

393-
394393
/////////////////////////////
395394
//
396395
// Set analytics
@@ -467,6 +466,23 @@ int user1_init_primitives(FTYPE (*prim)[NSTORE2][NSTORE3][NPR], FTYPE (*pstag)[N
467466

468467
if (bound_allprim(STAGEM1,0.0,prim,pstag,ucons, 1) >= 1) FAILSTATEMENT("init.c:init()", "bound_allprim()", 1);
469468

469+
470+
471+
/////////////////////////////
472+
//
473+
// Set EOSextras related to keeping table result consistent
474+
//
475+
///////////////////////////////
476+
477+
if(WHICHEOS==KAZFULL){
478+
FULLLOOP{
479+
// then store pressure
480+
// assume standard inversion at loc=CENT
481+
GLOBALMACP0A1(EOSextraglobal,i,j,k,PGASGLOBAL)=pressure_rho0_u(WHICHEOS,GLOBALMAC(EOSextraglobal,i,j,k),MACP0A1(prim,i,j,k,RHO),MACP0A1(prim,i,j,k,UU));
482+
}
483+
}
484+
485+
470486
// now fully bounded, initialize interpolations in case interpolate using prim/pstag data
471487
pre_interpolate_and_advance(prim);
472488

‎initboundcode/init.grb.c

+1-1
Original file line numberDiff line numberDiff line change
@@ -795,7 +795,7 @@ int init_global(void)
795795
gamideal=4.0/3.0; // ideal gas EOS index to use if off table
796796

797797
// below assumes tabulated EOS has minimum internal energy of zero, which is approximately correct currently.
798-
zerouuperbaryon=-TRUENUCLEAROFFSET; // definition of zero-point energy per baryon (i.e. such that internal energy cannot be lower than this times baryon density)
798+
zerouuperbaryon=-TRUENUCLEAROFFSETPRIMARY; // definition of zero-point energy per baryon (i.e. such that internal energy cannot be lower than this times baryon density)
799799

800800

801801
cooling=0;

‎jon_interp.c

+10-8
Original file line numberDiff line numberDiff line change
@@ -163,7 +163,7 @@ static void interp_readcommandlineargs(int argc, char *argv[])
163163
for(i=0;i<argc;i++){
164164
fprintf(stderr,"argv[%d]=%s\n",i,argv[i]);
165165
}
166-
fprintf(stderr,"args (argc=%d should be 27+ (26+ user args)): DATATYPE,INTERPTYPE,READHEADER,WRITEHEADER,oN1,oN2,oN3,refinefactor,filter,sigma,oldgridtype,newgridtype,nN1,nN2,nN3,startxc,endxc,startyc,endyc,startzc,endzc,Rin,Rout,R0,hslope,defcoord\n",argc) ;
166+
fprintf(stderr,"args (argc=%d should be 27+ (26+ user args)): DATATYPE,INTERPTYPE,READHEADER,WRITEHEADER,oN1,oN2,oN3,refinefactor,filter,sigma,oldgridtype,newgridtype,nN1,nN2,nN3,startxc,endxc,startyc,endyc,startzc,endzc,Rin,Rout,R0,hslope,defcoord,dofull2pi\n",argc) ;
167167

168168
fprintf(stderr,"DATATYPE:\n"
169169
"0=image (byte binary only, 1 column only)\n"
@@ -203,6 +203,7 @@ static void interp_readcommandlineargs(int argc, char *argv[])
203203
fprintf(stderr,"R0: Radial inner-grid enhancement factor\n");
204204
fprintf(stderr,"hslope: theta grid refinement factor\n");
205205
fprintf(stderr,"defcoord: which coordinate system (see coord.c)\n");
206+
fprintf(stderr,"dofull2pi: whether to do full 2pi or not (see coord.c)\n");
206207
fprintf(stderr,"extrapolate: 0 = no, 1 = yes\n");
207208
fprintf(stderr,"defaultvaluetype: 0 = min if scalar 0 if vector, 1 = min, 2 = max, 3 = 0.0, 4 = 1E35 for v5d missingdata\n");
208209
fprintf(stderr,"gdumpfilepathname : only if vector type\n");
@@ -323,6 +324,7 @@ static void interp_readcommandlineargs(int argc, char *argv[])
323324
sscanf(argv[i++],SCANARG,&R0) ;
324325
sscanf(argv[i++],SCANARG,&hslope) ;
325326
sscanf(argv[i++],"%d",&defcoord) ;
327+
sscanf(argv[i++],"%d",&dofull2pi) ;
326328

327329
// conditionally read-in things (all in or out)
328330
getgdump=0;
@@ -431,7 +433,7 @@ static void setup_zones(void)
431433
a=spin;
432434
Rhor=rhor_calc(0);
433435

434-
fprintf(stderr,"defcoord=%d Rin=%g R0=%g\n",defcoord,Rin,R0);
436+
fprintf(stderr,"defcoord=%d Rin=%g R0=%g dofull2pi=%d\n",defcoord,Rin,R0,dofull2pi);
435437

436438
fprintf(stderr,"startx,Xmax,dX: %g %g %g :: %g %g %g :: %g %g %g\n",startx[1],Xmax[1],dX[1],startx[2],Xmax[2],dX[2],startx[3],Xmax[3],dX[3]) ; fflush(stderr);
437439

@@ -627,7 +629,7 @@ static void input_header(void)
627629
// GODMARK
628630
// If using gammie.m's interpsingle, must keep interpsingle macro's header output up-to-date
629631
fscanf(stdin, SCANHEADER,
630-
&t,&totalsize[1],&totalsize[2],&totalsize[3],&startx[1],&startx[2],&startx[3],&dX[1],&dX[2],&dX[3],&readnstep,&gam,&spin,&R0,&Rin,&Rout,&hslope,&dt,&defcoord,&MBH,&QBH,&is,&ie,&js,&je,&ks,&ke);
632+
&t,&totalsize[1],&totalsize[2],&totalsize[3],&startx[1],&startx[2],&startx[3],&dX[1],&dX[2],&dX[3],&readnstep,&gam,&spin,&R0,&Rin,&Rout,&hslope,&dt,&defcoord,&MBH,&QBH,&is,&ie,&js,&je,&ks,&ke,&whichdump,&whichdumpversion,&numcolumns);
631633

632634

633635

@@ -641,7 +643,7 @@ static void input_header(void)
641643

642644
// print header from file
643645
fprintf(stderr, PRINTSCANHEADER,
644-
t,totalsize[1],totalsize[2],totalsize[3],startx[1],startx[2],startx[3],dX[1],dX[2],dX[3],readnstep,gam,spin,R0,Rin,Rout,hslope,dt,defcoord,MBH,QBH,is,ie,js,je,ks,ke);
646+
t,totalsize[1],totalsize[2],totalsize[3],startx[1],startx[2],startx[3],dX[1],dX[2],dX[3],readnstep,gam,spin,R0,Rin,Rout,hslope,dt,defcoord,MBH,QBH,is,ie,js,je,ks,ke,whichdump,whichdumpversion,numcolumns);
645647

646648

647649
}
@@ -658,8 +660,8 @@ static void output_header(void)
658660
//
659661
//////////////////////////
660662
fprintf(stderr,"header:\n");
661-
fprintf(stderr, "OLD: %22.16g :: %d %d %d :: %22.16g %22.16g %22.16g :: %22.16g %22.16g %22.16g :: %ld %22.16g %22.16g %22.16g %22.16g %22.16g %22.16g %22.16g %d %22.16g %22.16g %d %d %d %d %d %d\n",
662-
t,oN1,oN2,oN3,startx[1],startx[2],startx[3],dX[1],dX[2],dX[3],realnstep,gam,spin,R0,Rin,Rout,hslope,dt,defcoord,MBH,QBH,is,ie,js,je,ks,ke);
663+
fprintf(stderr, "OLD: %22.16g :: %d %d %d :: %22.16g %22.16g %22.16g :: %22.16g %22.16g %22.16g :: %ld %22.16g %22.16g %22.16g %22.16g %22.16g %22.16g %22.16g %d %22.16g %22.16g %d %d %d %d %d %d %d %d %d\n",
664+
t,oN1,oN2,oN3,startx[1],startx[2],startx[3],dX[1],dX[2],dX[3],realnstep,gam,spin,R0,Rin,Rout,hslope,dt,defcoord,MBH,QBH,is,ie,js,je,ks,ke,whichdump,whichdumpversion,numcolumns);
663665
fprintf(stderr, "NEW: %d %d %d :: %22.16g %22.16g %22.16g :: %22.16g %22.16g %22.16g\n",nN1,nN2,nN3,Xmax[1],Xmax[2],Xmax[3],fakedxc,fakedyc,fakedzc);
664666
fprintf(stderr, "OTHER: %22.16g %22.16g %22.16g %22.16g\n",fakeRin,dxc,dyc,dzc);
665667

@@ -668,8 +670,8 @@ static void output_header(void)
668670
if(DATATYPE==1){
669671
// print out a header
670672
ftemp=0.0;
671-
fprintf(stdout, "%22.16g %d %d %d %22.16g %22.16g %22.16g %22.16g %22.16g %22.16g %ld %22.16g %22.16g %22.16g %22.16g %22.16g %22.16g %22.16g %d %22.16g %22.16g %d %d %d %d %d %d\n",
672-
t, nN1, nN2, nN3, startxc, startyc, startzc, fakedxc,fakedyc,fakedzc,realnstep,gam,spin,ftemp,endxc,endyc,hslope,dt,defcoord,MBH,QBH,is,ie,js,je,ks,ke);
673+
fprintf(stdout, "%22.16g %d %d %d %22.16g %22.16g %22.16g %22.16g %22.16g %22.16g %ld %22.16g %22.16g %22.16g %22.16g %22.16g %22.16g %22.16g %d %22.16g %22.16g %d %d %d %d %d %d %d %d %d\n",
674+
t, nN1, nN2, nN3, startxc, startyc, startzc, fakedxc,fakedyc,fakedzc,realnstep,gam,spin,ftemp,endxc,endyc,hslope,dt,defcoord,MBH,QBH,is,ie,js,je,ks,ke,whichdump,whichdumpversion,numcolumns);
673675
}
674676
}
675677

0 commit comments

Comments
 (0)