Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add code for modeling a radiation pressure dominated disk #14

Open
wants to merge 14 commits into
base: master
Choose a base branch
from
Open

Add code for modeling a radiation pressure dominated disk #14

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
14 commits
Select commit Hold shift + click to select a range
e341ca5
Add code for a radiation pressure dominated disk
soumide1102 Mar 1, 2021
2b612f2
Fix a commented line
soumide1102 Mar 1, 2021
be323b4
Keep EOS_TYPE_GAMMA and break that into GAMMA_EOS = RADPRESS and GASP…
soumide1102 Mar 2, 2021
79aec3b
Fix the gamma option names consistent to be consistent
soumide1102 Mar 2, 2021
461f505
GAMMA_EOS --> EOS_GAMMA in all occurences
soumide1102 Mar 2, 2021
24aee2a
EOS_TYPE --> EOS everywhere
soumide1102 Mar 2, 2021
d136659
Update the remaining code with the two updated EOS_TYPE_GAMMA options
soumide1102 Mar 9, 2021
d51238c
Trying to fix ent uninitialized err from unittest
soumide1102 Mar 9, 2021
f42c723
Add a default value for EOS_GAMMA
soumide1102 Mar 9, 2021
a76212a
Multiplied in rad constant in rho eqn and added rad dominated option …
soumide1102 Mar 11, 2021
3390946
Forgot to add eos_gamma in the prev commit
soumide1102 Mar 11, 2021
06164ff
edits to make entropy be declared in problem.c and some debugging prints
soumide1102 Mar 18, 2021
87b77cb
switching back to correctly initializing entropy in problem.c instead…
soumide1102 May 4, 2021
6a33171
Calculate in CGS then divide by code units
soumide1102 May 25, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 9 additions & 3 deletions core/decs.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -201,7 +201,9 @@
#if EOS == EOS_TYPE_GAMMA
#define EOS_NUM_EXTRA (0)
#define POLYTROPE_FALLBACK (0)
#elif EPS == EOS_TYPE_POLYTROPE
#define GASPRESS (1)
#define RADPRESS (2)
#elif EOS == EOS_TYPE_POLYTROPE
Comment on lines +204 to +206
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes. This is what I was thinking. :)

#define EOS_NUM_EXTRA (0)
#define POLYTROPE_FALLBACK (0)
#elif EOS == EOS_TYPE_TABLE
Expand Down Expand Up @@ -417,6 +419,9 @@ extern double a;
#if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK
extern double gam;
#endif
//#if EOS == EOS_TYPE_TABLE || (EOS == EOS_TYPE_GAMMA && EOS_GAMMA == RADPRESS)
//extern double entropy;
//#endif
#if EOS == EOS_TYPE_POLYTROPE
extern double poly_K, poly_gam;
#endif
Expand All @@ -427,7 +432,7 @@ extern double M_unit;
extern double Reh;
extern double Risco;
#if NEED_UNITS
extern double mbh, Mbh, L_unit, T_unit, M_unit, RHO_unit, U_unit, B_unit;
extern double mbh, Mbh, L_unit, T_unit, M_unit, RHO_unit, U_unit, B_unit, TEMP_unit;
#endif
#if EOS == EOS_TYPE_TABLE
extern double TEMP_unit;
Expand Down Expand Up @@ -787,7 +792,8 @@ double EOS_Theta_unit();
#if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK
double EOS_Gamma_pressure_rho0_u(double rho, double u);
double EOS_Gamma_pressure_rho0_w(double rho, double w);
double EOS_Gamma_entropy_rho0_u(double rho, double u);
double EOS_Gamma_Gaspress_entropy_rho0_u(double rho, double u);
double EOS_Gamma_Radpress_entropy_rho0_u(double rho, double u);
double EOS_Gamma_enthalpy_rho0_u(double rho, double u);
double EOS_Gamma_adiabatic_constant_rho0_u(double rho, double u);
double EOS_Gamma_sound_speed_rho0_u(double rho, double u);
Expand Down
5 changes: 4 additions & 1 deletion core/defs.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -101,6 +101,9 @@ char opac_file[STRLEN];
#if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK
double gam;
#endif
//#if EOS == EOS_TYPE_TABLE || (EOS == EOS_TYPE_GAMMA && EOS_GAMMA == RADPRESS)
//double entropy;
//#endif
#if EOS == EOS_TYPE_POLYTROPE
double poly_K, poly_gam;
#endif
Expand All @@ -114,7 +117,7 @@ double Reh;
double Risco;

#if NEED_UNITS
double mbh, Mbh, L_unit, T_unit, M_unit, RHO_unit, U_unit, B_unit;
double mbh, Mbh, L_unit, T_unit, M_unit, RHO_unit, U_unit, B_unit, TEMP_unit;
#endif

#if EOS == EOS_TYPE_TABLE
Expand Down
6 changes: 5 additions & 1 deletion core/eos.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -160,7 +160,11 @@ double EOS_enthalpy_rho0_u(double rho, double u, const double *extra) {
double EOS_entropy_rho0_u(double rho, double u, const double *extra) {
double ent;
#if EOS == EOS_TYPE_GAMMA
ent = EOS_Gamma_entropy_rho0_u(rho, u);
#if EOS_GAMMA == GASPRESS
ent = EOS_Gamma_Gaspress_entropy_rho0_u(rho, u);
#elif EOS_GAMMA == RADPRESS
ent = EOS_Gamma_Radpress_entropy_rho0_u(rho, u);
#endif
#elif EOS == EOS_TYPE_POLYTROPE
ent = EOS_Poly_entropy_rho0_u(rho, u, poly_K, poly_gam);
#elif EOS == EOS_TYPE_TABLE
Expand Down
23 changes: 20 additions & 3 deletions core/eos_gamma.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,10 +17,19 @@ double EOS_Gamma_pressure_rho0_w(double rho, double w) {
return ((w - rho) * (gam - 1.) / gam);
}

double EOS_Gamma_entropy_rho0_u(double rho, double u) {
double EOS_Gamma_Gaspress_entropy_rho0_u(double rho, double u) {
return (gam - 1.0) * EOS_Gamma_adiabatic_constant_rho0_u(rho, u);
}

double EOS_Gamma_Radpress_entropy_rho0_u(double rho, double u) {
//double a = AR * pow(T_unit,2) * L_unit * pow(TEMP_unit, 4) * pow(M_unit, -2);
double u_cgs = u * U_unit;
double rho_cgs = rho * RHO_unit;
double entropy_cgs = pow((64./3) * AR * (pow(u_cgs, 3)*pow(gam, 3)/pow(rho_cgs, 4)) * pow((gam -1)/gam, 3), 1./4);
//return pow((64./3) * a * (pow(u, 3)*pow(gam, 3)/pow(rho, 4)) * pow((gam -1)/gam, 3), 1./4);
return entropy_cgs / (KBOL / MP);
}

double EOS_Gamma_enthalpy_rho0_u(double rho, double u) { return rho + u * gam; }

double EOS_Gamma_adiabatic_constant_rho0_u(double rho, double u) {
Expand Down Expand Up @@ -62,8 +71,16 @@ double EOS_Gamma_sound_speed_rho0_u(double rho, double u) {
double EOS_Gamma_u_press(double press) { return press / (gam - 1.); }

double EOS_Gamma_temp(double rho, double u) {
// return TEMP_unit*(gam - 1.)*u/rho;
return (gam - 1.) * u / rho;
#if EOS_GAMMA == RADPRESS
double press_cgs = (gam - 1.) * u * U_unit;
double temp_cgs = pow((3 * press_cgs)/AR, 1./4);
double temp = temp_cgs * KBOL / MEV;
//double temp = pow((3 * press), 1./4);
#else
// return TEMP_unit*(gam - 1.)*u/rho;
double temp = (gam - 1.) * u / rho;
#endif
return temp;
}

#if RADIATION
Expand Down
8 changes: 6 additions & 2 deletions core/fixup.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -229,9 +229,9 @@ void fixup1zone(
fixup_passive(i, j, k, pv, pv_prefloor);
#endif

#if ELECTRONS && EOS == EOS_TYPE_GAMMA
#if ELECTRONS && EOS == EOS_TYPE_GAMMA && GAMMA_EOS == GASPRESS
// Reset entropy after floors
pv[KTOT] = EOS_Gamma_entropy_rho0_u(pv[RHO], pv[UU]);
pv[KTOT] = EOS_Gamma_Gaspress_entropy_rho0_u(pv[RHO], pv[UU]);

// Set KTOTMAX to 3 by controlling u, to avoid anomalous cooling from funnel
// wall
Expand All @@ -240,6 +240,10 @@ void fixup1zone(
pv[UU] = KTOTMAX * pow(pv[RHO], gam) / (gam - 1.);
pv[KTOT] = KTOTMAX;
}

#elif ELECTRONS && EOS == EOS_TYPE_GAMMA && GAMMA_EOS == RADPRESS
// Reset entropy after floors
pv[KTOT] = EOS_Gamma_Radpress_entropy_rho0_u(pv[RHO], pv[UU]);
#endif // ELECTRONS

// Limit gamma with respect to normal observer
Expand Down
4 changes: 4 additions & 0 deletions core/input.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -95,6 +95,10 @@ void set_core_params() {
#if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK
set_param("gam", &gam);
#endif
//#if EOS == EOS_TYPE_TABLE || (EOS == EOS_TYPE_GAMMA && EOS_GAMMA == RADPRESS)
// set_param("entropy", &entropy);
//#endif

#if EOS == EOS_TYPE_GAMMA
sprintf(eos, "GAMMA");
#elif EOS == EOS_TYPE_POLYTROPE
Expand Down
8 changes: 5 additions & 3 deletions core/io.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -730,7 +730,7 @@ void dump() {
#endif // RADIATION
#endif // METRIC

#if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK
#if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK
WRITE_HDR(gam, TYPE_DBL);
#endif // GAMMA
#if EOS == EOS_TYPE_POLYTROPE
Expand All @@ -753,6 +753,7 @@ void dump() {
WRITE_HDR(T_unit, TYPE_DBL);
WRITE_HDR(U_unit, TYPE_DBL);
WRITE_HDR(B_unit, TYPE_DBL);
WRITE_HDR(TEMP_unit, TYPE_DBL);
#if EOS == EOS_TYPE_TABLE
WRITE_HDR(TEMP_unit, TYPE_DBL);
#endif // EOS_TYPE_TABLE
Expand Down Expand Up @@ -1099,7 +1100,7 @@ void restart_write(int restart_type) {
#endif // METRIC

// EOS
#if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK
#if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK
WRITE_HDR(gam, TYPE_DBL);
#endif // GAMMA EOS
#if EOS == EOS_TYPE_POLYTROPE
Expand All @@ -1117,6 +1118,7 @@ void restart_write(int restart_type) {
WRITE_HDR(T_unit, TYPE_DBL);
WRITE_HDR(U_unit, TYPE_DBL);
WRITE_HDR(B_unit, TYPE_DBL);
WRITE_HDR(TEMP_unit, TYPE_DBL);
#if EOS == EOS_TYPE_TABLE
WRITE_HDR(TEMP_unit, TYPE_DBL);
#endif // EOS_TYPE_TABLE
Expand Down Expand Up @@ -1381,7 +1383,7 @@ void restart_read(char *fname) {

READ_HDR(cour, TYPE_DBL);

#if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK
#if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK
READ_HDR(gam, TYPE_DBL);
#endif // EOS

Expand Down
5 changes: 5 additions & 0 deletions core/main.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -59,6 +59,11 @@ int main(int argc, char *argv[]) {
#endif
#if EOS == EOS_TYPE_GAMMA
fprintf(stdout, " * IDEAL GAS EOS *\n");
#if EOS_GAMMA == GASPRESS
fprintf(stdout, " * GAS-PRESSURE DOMINATED IDEAL GAS EOS *\n");
#elif EOS_GAMMA == RADPRESS
fprintf(stdout, " * RADIATION-PRESSURE DOMINATED IDEAL GAS EOS *\n");
#endif
#endif
#if RADIATION == RADTYPE_LIGHT
fprintf(stdout, " * RADIATION: PHOTON TRANSPORT *\n");
Expand Down
2 changes: 1 addition & 1 deletion core/phys.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -211,7 +211,7 @@ void mhd_vchar(double *Pr, struct of_state *q, struct of_geom *geom, int js,
bsq = dot(q->bcon, q->bcov);
rho = fabs(Pr[RHO] + SMALL);
u = Pr[UU];
#if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK
#if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK
u = fabs(u + SMALL);
#endif
ef = EOS_enthalpy_rho0_u(rho, u, extra);
Expand Down
3 changes: 3 additions & 0 deletions core/util.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -145,6 +145,9 @@ void set_units() {
B_unit = CL * sqrt(4. * M_PI * RHO_unit);
#if EOS == EOS_TYPE_TABLE
TEMP_unit = MEV;
#elif EOS == EOS_TYPE_GAMMA
//TEMP_unit = GNEWT * KBOL / pow(CL, 4);
TEMP_unit = (1. / KBOL) * (M_unit * pow(L_unit, 2.) / pow(T_unit, 2.));
#endif

#if RADIATION
Expand Down
2 changes: 1 addition & 1 deletion core/utop.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -161,7 +161,7 @@ int Utoprim(double U[NVAR], struct of_geom *geom, double prim[NVAR]) {
// Return without updating non-B primitives
if (rho0 < 0)
return (5);
#if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK
#if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK
if (u < 0)
return (5);
#endif
Expand Down
50 changes: 45 additions & 5 deletions prob/torus_cbc/build.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,8 @@
PROB = 'torus_cbc'

DO_GAMMA = '-gamma' in sys.argv # No table
GAMMA_GASPRESS = '-gammagaspress' in sys.argv # No table
GAMMA_RADPRESS = '-gammaradpress' in sys.argv # No table
GAMTABLE = '-gamtable' in sys.argv # fake table
RELTABLE = '-reltable' in sys.argv or (not DO_GAMMA or GAMTABLE)
INITIAL = '-initial' in sys.argv # initial data only
Expand Down Expand Up @@ -51,6 +53,7 @@
KILL = '-kill' in sys.argv # kill all packets

ENT_FROM_CLI = '-ent' in sys.argv
KAPPA_FROM_CLI = '-kappa' in sys.argv
RHO_FROM_CLI = '-rho' in sys.argv
RIN_FROM_CLI = '-rin' in sys.argv
RMAX_FROM_CLI = '-rmax' in sys.argv
Expand Down Expand Up @@ -83,6 +86,9 @@
if any([N1TOT_FROM_CLI, N2TOT_FROM_CLI, N3TOT_FROM_CLI]) and not N1N2N3TOT_FROM_CLI:
raise ValueError("Must turn on the -n1n2n3tot flag if inputting -n1tot, -n2tot, -n3tot")

if any([GAMMA_GASPRESS, GAMMA_RADPRESS]) and not DO_GAMMA:
raise ValueError("Must turn on the -gamma flag if inputting -gammagaspress or -gammaradpress")

if NOB:
BFIELD = "none"
elif TOROIDALB:
Expand Down Expand Up @@ -146,12 +152,43 @@
if USE_TABLE:
bhl.report_var('Disk Ye',YE)

if USE_TABLE:
EOS_TYPE = "EOS_TYPE_TABLE"
if DO_GAMMA:
EOS_TYPE = "EOS_TYPE_GAMMA"
if GAMMA_RADPRESS:
EOS_GAMMA = "RADPRESS"
GAMMA = 4./3.
elif GAMMA_GASPRESS:
EOS_GAMMA = "GASPRESS"
GAMMA = 5./3.
else:
print('Using the default gas-pressure dominated ideal gas EoS.',
'You can choose gas-pressure dominated (add flag -gammagaspress)',
'radiation-pressure dominated (add flag -gammaradpress).')
EOS_GAMMA = "GASPRESS"
GAMMA = 4./3.
else:
raise ValueError("Bad EOS chosen")

if ENT_FROM_CLI:
ENTROPY = float(sys.argv[sys.argv.index('-ent') + 1])
else:
ENTROPY = 4
if EOS_TYPE == "EOS_TYPE_GAMMA" and EOS_GAMMA == "RADPRESS":
print('Using the default initial entropy = 4.',
'Entropy is a required parameter in setting up radiation-pressure dominated disks.')

bhl.report_var('ENTROPY',ENTROPY)

if KAPPA_FROM_CLI:
KAPPA_EOS = float(sys.argv[sys.argv.index('-kappa') + 1])
else:
KAPPA_EOS = 1.e-3
if EOS_TYPE == "EOS_TYPE_GAMMA" and EOS_GAMMA == "GASPRESS":
print('Using the default initial kappa = 1.e-3.',
'kappa is a required parameter in setting up gas-pressure dominated disks.')

if CLASSIC: # classic harm disk
# Rmax and rho fixed
bhl.report_var('DISK_TYPE','CLASSIC')
Expand Down Expand Up @@ -211,8 +248,6 @@
# use selection instead
Rout_rad = ENTROPY*ceil(Rmax) # not safe to use 3x
Rout_vis = Rout
GAMMA = 13./9.
KAPPA = 1.e-3
L_unit = cgs['GNEWT']*cgs['MSOLAR']*MBH/(cgs['CL']**2)
M_UNIT = RHO_unit*(L_unit**3)

Expand Down Expand Up @@ -314,7 +349,6 @@
if not NEUTRINOS:
NPH_PER_PROC = 0.0

EOS = "EOS_TYPE_TABLE" if USE_TABLE else "EOS_TYPE_GAMMA"
NVP0 = 3
NVAR_PASSIVE = NVP0 if USE_TABLE else 0
#GAMMA_FALLBACK = True
Expand Down Expand Up @@ -389,7 +423,10 @@
bhl.config.set_cparm('QUADRANT_SYMMETRY', QUAD)

# EOS
bhl.config.set_cparm("EOS", EOS)
bhl.config.set_cparm("EOS", EOS_TYPE)
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So one small change I was thinking of: Here the name of the parameter is EOS which saves which type of EoS is going to be used. Also, around line 100, there is

if RELTABLE:
    bhl.report_var('EOS','RELTABLE')
elif GAMTABLE:
    bhl.report_var('EOS','GAMTABLE')
else:
    bhl.report_var('EOS','GAMMA')

where the parameter name is EOS. I was wondering whether the latter parameter can be changed to something like EOS_VAL to keep the two parameter names different?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

no those should be the same. report_var doesn't do anything internal in the code. It prints a choice of variable to the build output.

bhl.config.set_cparm('OUTPUT_EOSVARS', True)
if EOS_TYPE == 'EOS_TYPE_GAMMA':
bhl.config.set_cparm("EOS_GAMMA", EOS_GAMMA)
bhl.config.set_cparm('NVAR_PASSIVE', NVAR_PASSIVE)
bhl.config.set_cparm('GAMMA_FALLBACK', GAMMA_FALLBACK)

Expand Down Expand Up @@ -466,7 +503,10 @@
bhl.config.set_rparm('lrho_guess','double',LRHO_GUESS)
if USE_GAMMA or GAMMA_FALLBACK:
bhl.config.set_rparm('gam', 'double', default = GAMMA)
bhl.config.set_rparm("kappa", "double", default = KAPPA)
if EOS_TYPE == "EOS_TYPE_GAMMA" and EOS_GAMMA == "RADPRESS":
bhl.config.set_rparm('entropy','double',ENTROPY)
elif EOS_TYPE == "EOS_TYPE_GAMMA" and EOS_GAMMA == "GASPRESS":
bhl.config.set_rparm("kappa_eos", "double", default = KAPPA_EOS)

# Opacities
if FORTRAN or HDF:
Expand Down
Loading