Cholla unit system

[brantr]

One issue that I need to address in resolving problems with the cosmology runs and in merging in the RT is the unit system. This issue could affect others, so I thought I'd describe what I'm thinking and request discussion.

There are several of the unit systems defined in the code. Many of these units are defined in global.h. However, there are internal inconsistencies. For instance:

MYR                                  3.153600000e+13. [myr in s]
TIME_UNIT                            3.155690000e+10. [kyr in s]

Neither of these are what I expected, although TIME_UNIT is close.

Also, this:

#define MP          1.672622e-24  // mass of proton, grams
#define MH       1.67262171e-24    // Mass of hydrogen [g]

Also GN is defined in kpc^3 / M_sun / kyr^2 but there's a commented-out CGS version.

There are other unit system questions I have with the cosmology and chemistry unit systems. Some units are set inline in functions.

Here are the things I would like to pursue:

1. Refining the unit systems defined in global.h. Clearly this could impact everyone and unit tests, so this would not be taken lightly and would need input from everyone. Relabel all units that are in CGS with a _CGS suffix. Label everything unit in the kpc|km|s|Msun system and the kpc|kyr|Msun system with distinct suffixes.

2. Writing a Grid3D member function that will print the actual unit system to stdout at the beginning of a simulation, called in main(). This could be couched in a preproc def and set at compile time to print or not, but really it would be helpful. Or it could be written to a separate text output file.

3. Write functions for cosmology and chemistry where all the units for those modules are set in function members of Cosmology and Chem.

4. Write functions that print the cosmology and chemistry units when they are defined, which would be called by the grid member function.

While I listed 1-4, 2-4 is actually a tractable and small change that should impact no one. I am planning to submit a PR that addresses 2-4, but 1) probably requires discussion and so I thought I'd submit this issue. Thanks!

[evaneschneider]

Thanks for starting this discussion, Brant. I agree that the current unit system needs reform. Your approach for 2 - 4 seems very reasonable. Regarding 1, I think a lot of the duplicate units were added when the cosmology / chemistry code was added, so we will probably need to update test data, initial conditions, the older cooling code, etc. if the plan is to reconcile to the newer units. I am happy to help with that. The only units that should definitely stay defined are the ones labeled things like TIME_UNIT, LENTH_UNIT, etc. These have always been intended to be the internal units that Cholla uses to convert between code units and CGS, and they are currently output in the hdf5 header files and used by other analysis code (for example, the Cholla yt frontend). That said, if the actual values are off (I have no idea why kyr is defined that way instead of as 3.1536e10, for example), we should correct them. (And again, I'm happy to help with that, as needed.)

[brantr]

Thanks @evaneschneider .

I will work on the functions for items 2-4 and submit a PR for them.

Regarding TIME_UNIT and kyr, it's just a question of whether to pick 365 or ~365.24... or 365.25 as a "year". Note that NIST defines a lightyear with 1 day = 86400s and a Juilan century as 36525 days, which would be 3.15576e10 seconds in a kiloyear. This definition is consistent with the NIST value for a light year of 9.46073e17 cm when defining c = 29979245800 cm/s.

https://www.nist.gov/pml/special-publication-811/nist-guide-si-appendix-b-conversion-factors/nist-guide-si-appendix-b8#L

[mabruzzo]

There are other unit system questions I have with the cosmology and chemistry unit systems. Some units are set inline in functions.

Interjecting here: does this mean we track separate unit-systems for chemistry and cosmology? If so, does that mean that Cholla's cosmological unit-systems is incompatible with Grackle's cosmological unit system? (Figuring out whether Grackle's unit-system with any cosmological unit-systems other than Enzo's cosmological unit system is something I would like to figure out -- grackle-project/grackle#216 ).

Write functions for cosmology and chemistry where all the units for those modules are set in function members of Cosmology and Chem.

This probably isn't an issue for right now, but would it actually make more sense to unify the Cosmology and Hydro unit systems? From a performance standpoint, there isn't a lot of advantage to tracking TIME_UNIT, LENTH_UNIT, etc. as macros. If we instead created an object that tracked the simulation's initial units and then provided a function that could compute the units at an arbitrary time (in non-cosmological sims this wouldn't do any work), that would probably bring us a lot closer to not needing to separately recompile Cholla for cosmological/non-cosmological problems. (If you guys don't think this is valuable, that's totally fine with me -- I realize my background with Enzo-E may skew my perspective here).

[bvillasen]

Grackle uses GCS. You need to tell geackle how to convert from your units to GCS.
From Grid3D::Do_Cooling_Step_Grackle():

Real kpc_cgs = KPC_CGS;
// Update the units conversion
Cool.units.a_value       = Cosmo.current_a / Cool.units.a_units;
Cool.units.density_units = Cool.dens_to_CGS / Cosmo.current_a / Cosmo.current_a / Cosmo.current_a;
Cool.units.length_units  = kpc_cgs / Cosmo.cosmo_h * Cosmo.current_a;

Copy_Fields_To_Grackle();

Real dt_cool = Cosmo.dt_secs;
chprintf(" dt_cool: %e s\n", dt_cool);
if (solve_chemistry(&Cool.units, &Cool.fields, dt_cool / Cool.units.time_units) == 0) {
  chprintf("GRACKLE: Error in solve_chemistry.\n");
  return;
}

You see that we need to convert from comoving to physical units.
Also, I remember that gracke takes the specific internal energy U/rho while we track U (when using dual energy).
So yes, there are unit conversion needed to use grackle. Or at least that was the case ~4 years ago, not sure if it's still the case.
Happy to talk more if you want, just note that I might not remember much :)

[mabruzzo]

So, as of March, I'm now co-manager of Grackle (alongside Britton Smith) :). Making Grackle more user-friendly has been one of my goals.

Grackle uses GCS. You need to tell geackle how to convert from your units to GCS.
...
You see that we need to convert from comoving to physical units.

Right, so Grackle actually supports 2 types of unit-systems.

1. It supports a proper unit-system (where you specify the conversion between code_units and CGS).
   • It sounds like Cholla currently uses Grackle with this approach.
   • As you say, cosmological Cholla sims needs to converts the units of all Cholla fields to and from proper units in this scenario.
2. It also supports a comoving coordinate system.
   • As you may recall, during initialization, Grackle pre-computes certain cooling rates and stores these values internally. Essentially, Grackle requires that the comoving coordinate system is defined in such a way that the conversion between the code_units for each of these internal quantities and cgs units has no dependence on redshift
   • If a cosmological simulation code (like Enzo-Classic or Enzo-E) employs a comoving unit-system the satisfies this invariant, then you don't need to convert the units of the simulation's fields to/from comoving units before/after calling Grackle.

The reason I bring this all up: I want to make Grackle easier to use (and harder to misuse).

• It is currently extremely easy to misuse Grackle with comoving unit support.
• In fact, I'm fairly confident that Grackle's comoving unit support is only compatible with a comoving unit-system just like Enzo's.
• If that is true, then I have proposed a (backwards-compatible) modification to Grackle's API to make it harder to misuse Grackle (we might eventually be able build in support for other comoving unit-systems).
• The problem is that Britton and I aren't absolutely certain Grackle can't currently be used with a different comoving system. From this conversation, it doesn't sound like Grackle's comoving support is compatible with the Ramses/Cholla comoving unit-system (or it hasn't been rigorously tested),

Also, I remember that gracke takes the specific internal energy U/rho while we track U (when using dual energy).

Yep! That's exactly right! This is very much a product of the fact that Grackle was extracted from Enzo-Classic.

So yes, there are unit conversion needed to use grackle. Or at least that was the case ~4 years ago, not sure if it's still the case.

Things definitely haven't really changed in this area over the last 4 years.

Happy to talk more if you want, just note that I might not remember much :)

Chatting more could be useful. There's also ongoing efforts to get Grackle onto GPUs (which obviously involves porting a lot of Fortran to C/C++). We have looked at some of the stuff that you have done. But, we're trying to walk a tricky line of maintaining entirely consistency with existing behavior, maintaining support for simulation-codes on CPUs (without needing to maintain 2 copies of the codebase), and ideally we want to maintain the existing API.

[bvillasen]

Nice, I didn't know you had such an intimate relationship with Grackle.
Yes, I remember that figuring out the correct units to pass to Grackle took much longer than adding Grackle into Cholla, so I agree that improving the unit system will make it much more user-friendly. I remember that I was confused because I have units.comoving_coordinates = 1; when I initialize grackle, but I still had to convert to physical units for the other unit conversions. I probably tried many combinations and that is the one that gave the correct result so I just sticked to that. Anyways, I'm glad you are working on improving this.

I'm not sure if you are aware, but for us, moving the H+He+UVB stuff that grackle does to the GPU gave a huge speedup, I want to say over 100x, and this was on Summit, on Frontier is probably more since there are fewer CPU cores. Before moving the Grackle stuff, over 60% of our run time on summit was spent on Grackle calculations, after moving to the GPU, the chemistry became a very small fraction of the run time.

And honestly I don't think that porting Grackle to GPUs would be that complicated either (at least not the primordial chemistry stuff, not sure about the more sophisticated stuff). The CPU arrays and GPU arrays can be maintained equally, and you could put the whole grackle computation for each cell in a single host device function and just call in in a kernel launch if you are running on the the GPU or in a loop if you are running on the CPU. That way you just maintain the content of that function and it should work for both... It's possible that I do not fully remember the whole complexity, but I can tell you from experience that it wasn't that difficult to reproduce that grackle functionality on the GPU.

If you think it would be useful, I'm happy to chat more about it, and maybe I can give some suggestions; just ping me on Slack.
And yes, please get rid of that Fortran :)