Add code for modeling a radiation pressure dominated disk #14

soumide1102 · 2021-03-01T20:56:44Z

This PR adds code for setting up a radiation pressure dominated disk.

For a gamma disk, the code now can either set up a radiation pressure dominated or a gas pressure dominated case. For doing that the user can input -gammagaspress or gammaradpress. -gammagaspress will follow the old code path for gamma. -gammaradpress will follow the radiation pressure code path added in this PR.
For the radiation-pressure dominated case, rho and u in problem.c is calculated using an inputted entropy.
Also modified the other files in nublight to incorporate the radiation pressure dominated condition.
The kappa parameter saved in build.py is now corrected to kappa_eos

In calculating rho for the radiation pressure dominated case, there is a multiplication factor a that is currently not included in problem.c . Should hdf5_to_dict be the place to factor that in when loading the results and converting to the desired units?

soumide1102 · 2021-03-01T23:00:21Z

core/fixup.c

+#if ELECTRONS && EOS == EOS_TYPE_GAMMA_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]);


I also wanted to check with you what this piece of code is doing. In the current modifications, I changed the name of the functions to the updated names specifying gas pressure. If this is going into evolving the simulation, should there also be another if statement for the radiation pressure dominated case as well?

This does go into evolving the simulation, yes. But this code path will never be reached. The ELECTRONS machinery is backwards compatibility with ebhlight's two-temperature gas for AGN disks. I would add another if case, but it's never going to come up.

With EOS_TYPE_GAMMA having two sub-cases, this can go to EOS_entropy_rho0_u or something

Yurlungur

Overall I like the strategy here... But something seems off as far as the code paths. You only replaced the entropy equation for the ideal gas equation of state, but you're treating it as an entirely different equation of state everywhere in the code. I think this is only working because you set all 3 macros to resolve to the same integer.

I would suggest a new structure. Something like:

The ideal gas eqation of state is triggered for EOS == EOS_TYPE_GAMMA
The gamma EOS has two new options GAMMA_GASPRESS and GAMMA_RADPRESS or something.
Those new macros are set as a new cparm in the build fie. We probably also need to set a new default value in script/configure.py

Yurlungur · 2021-03-02T01:02:11Z

core/decs.h

+#define EOS_TYPE_GAMMA_GASPRESS (0)
+#define EOS_TYPE_GAMMA_RADPRESS (0)


These all can't be zero. That means EOS_TYPE_GAMMA_GASPRESS and EOS_TYPEGAMMA_RADPRESS will resolve to the same number. That's probably not your intent. Each EOS type needs to resolve to a different number in the macro.

Yurlungur · 2021-03-02T01:06:41Z

core/main.c

@@ -60,6 +60,12 @@ int main(int argc, char *argv[]) {
    #if EOS == EOS_TYPE_GAMMA
    fprintf(stdout, "          *  IDEAL GAS EOS                                           *\n");
    #endif
+    #if EOS == EOS_TYPE_GAMMA_GASPRESS
+    fprintf(stdout, "          *  GAS-PRESSURE DOMINATED EOS                    *\n");


GAS-PRESSURE DOMINATED IDEAL GAS EOS

Yurlungur · 2021-03-02T01:06:52Z

core/main.c

+    fprintf(stdout, "          *  GAS-PRESSURE DOMINATED EOS                    *\n");
+    #endif
+    #if EOS == EOS_TYPE_GAMMA_RADPRESS
+    fprintf(stdout, "          *  RADIATION-PRESSURE DOMINATED EOS              *\n");


RADIATION-PRESSURE DOMINATED IDEAL GAS EOS

Yurlungur · 2021-03-02T01:13:46Z

prob/torus_cbc/build.py

@@ -211,8 +232,7 @@
 # use selection instead
 Rout_rad = ENTROPY*ceil(Rmax) # not safe to use 3x
 Rout_vis = Rout
-GAMMA = 13./9.
-KAPPA = 1.e-3
+GAMMA = 4./3.


This should be if-cased on whether you're radiation or gas pressure dominated. 4./3. is appropriate for radiation pressure, not gas pressure.

Yurlungur · 2021-03-02T01:14:40Z

prob/torus_cbc/build.py

+elif DO_GAMMA_RADPRESS:
+    EOS_TYPE = "EOS_TYPE_GAMMA_RADPRESS"
+elif DO_GAMMA_GASPRESS:
+    EOS_TYPE = "EOS_TYPE_GAMMA_GASPRESS"


now there should be a throw if it's not one a valid EOS type. Let's catch that in case someone somehow doesn't set one.

soumide1102 · 2021-03-02T03:40:54Z

@Yurlungur Does this look like what you were suggesting? I will then update the other scripts accordingly.

soumide1102 · 2021-03-02T03:57:35Z

prob/torus_cbc/build.py

@@ -408,7 +423,9 @@
 bhl.config.set_cparm('QUADRANT_SYMMETRY', QUAD)

 # EOS
-bhl.config.set_cparm("EOS_TYPE", EOS_TYPE)
+bhl.config.set_cparm("EOS", EOS_TYPE)


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?

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.

Yurlungur · 2021-03-02T16:25:03Z

core/decs.h

+#define GASPRESS (1)
+#define RADPRESS (2)
+#elif EOS == EOS_TYPE_POLYTROPE


Yes. This is what I was thinking. :)

Yurlungur · 2021-03-02T16:25:23Z

core/eos.c

+#if EOS == EOS_TYPE_GAMMA_GASPRESS
+  ent = EOS_Gamma_Gaspress_entropy_rho0_u(rho, u);
+#elif EOS == EOS_TYPE_GAMMA_RADPRESS
+  ent = EOS_Gamma_Radpress_entropy_rho0_u(rho, u);


Should now be sub-cases of EOS_TYPEGAMMA

Yurlungur · 2021-03-02T16:26:26Z

core/fixup.c

+#if ELECTRONS && EOS == EOS_TYPE_GAMMA_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]);


With EOS_TYPE_GAMMA having two sub-cases, this can go to EOS_entropy_rho0_u or something

Yurlungur · 2021-03-02T16:27:04Z

core/io.c

@@ -730,7 +730,7 @@ void dump() {
 #endif // RADIATION
 #endif // METRIC

-#if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK
+#if EOS == EOS_TYPE_GAMMA || EOS_TYPE_GAMMA_GASPRESS || EOS_TYPE_GAMMA_RADPRESS || GAMMA_FALLBACK  


can go back to just EOS == EOS_TYPE_GAMMA

Yurlungur · 2021-03-02T16:27:09Z

core/io.c

@@ -1099,7 +1099,7 @@ void restart_write(int restart_type) {
 #endif // METRIC

 // EOS
-#if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK
+#if EOS == EOS_TYPE_GAMMA || EOS_TYPE_GAMMA_GASPRESS || EOS_TYPE_GAMMA_RADPRESS || GAMMA_FALLBACK  


can go back to just EOS == EOS_TYPE_GAMMA

Yurlungur · 2021-03-02T16:27:13Z

core/io.c

@@ -1381,7 +1381,7 @@ void restart_read(char *fname) {

  READ_HDR(cour, TYPE_DBL);

-#if EOS == EOS_TYPE_GAMMA || GAMMA_FALLBACK
+#if EOS == EOS_TYPE_GAMMA || EOS_TYPE_GAMMA_GASPRESS || EOS_TYPE_GAMMA_RADPRESS || GAMMA_FALLBACK  


can go back to just EOS == EOS_TYPE_GAMMA

Yurlungur · 2021-03-09T18:54:43Z

@soumide1102 this is looking like it's moving along... although I see tests are failing. How close are you to being ready for another review?

…RESS

…in EOS_Gamma_temp

… of globally + experimenting with entropy CU

soumide1102 requested a review from jonahm-LANL March 1, 2021 20:56

soumide1102 commented Mar 1, 2021

View reviewed changes

Yurlungur requested changes Mar 2, 2021

View reviewed changes

soumide1102 commented Mar 2, 2021

View reviewed changes

Yurlungur reviewed Mar 2, 2021

View reviewed changes

soumide1102 added 7 commits March 9, 2021 14:10

Add code for a radiation pressure dominated disk

e341ca5

Fix a commented line

2b612f2

Keep EOS_TYPE_GAMMA and break that into GAMMA_EOS = RADPRESS and GASP…

be323b4

…RESS

Fix the gamma option names consistent to be consistent

79aec3b

GAMMA_EOS --> EOS_GAMMA in all occurences

461f505

EOS_TYPE --> EOS everywhere

24aee2a

Update the remaining code with the two updated EOS_TYPE_GAMMA options

d136659

soumide1102 force-pushed the rad_press_torus1 branch from 860cd19 to d136659 Compare March 9, 2021 21:18

soumide1102 added 7 commits March 9, 2021 14:26

Trying to fix ent uninitialized err from unittest

d51238c

Add a default value for EOS_GAMMA

f42c723

Multiplied in rad constant in rho eqn and added rad dominated option …

a76212a

…in EOS_Gamma_temp

Forgot to add eos_gamma in the prev commit

3390946

edits to make entropy be declared in problem.c and some debugging prints

06164ff

switching back to correctly initializing entropy in problem.c instead…

87b77cb

… of globally + experimenting with entropy CU

Calculate in CGS then divide by code units

6a33171

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add code for modeling a radiation pressure dominated disk #14

Add code for modeling a radiation pressure dominated disk #14

soumide1102 commented Mar 1, 2021

soumide1102 Mar 1, 2021

Yurlungur Mar 2, 2021

Yurlungur Mar 2, 2021

Yurlungur left a comment

Yurlungur Mar 2, 2021

Yurlungur Mar 2, 2021

Yurlungur Mar 2, 2021

Yurlungur Mar 2, 2021

Yurlungur Mar 2, 2021

soumide1102 commented Mar 2, 2021

soumide1102 Mar 2, 2021

Yurlungur Mar 9, 2021

Yurlungur Mar 2, 2021

Yurlungur Mar 2, 2021

Yurlungur Mar 2, 2021

Yurlungur Mar 2, 2021

Yurlungur Mar 2, 2021

Yurlungur Mar 2, 2021

Yurlungur commented Mar 9, 2021

		#define EOS_TYPE_GAMMA_GASPRESS (0)
		#define EOS_TYPE_GAMMA_RADPRESS (0)

Add code for modeling a radiation pressure dominated disk #14

Are you sure you want to change the base?

Add code for modeling a radiation pressure dominated disk #14

Conversation

soumide1102 commented Mar 1, 2021

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Yurlungur left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

soumide1102 commented Mar 2, 2021

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Yurlungur commented Mar 9, 2021