Auto-computing the Required Units

[mabruzzo]

In my opinion, one of the most difficult things about integrating Grackle into a simulation-code is its handling of the unit-system.

I think part of what makes it so confusing is that we place all the responsibility on the user for making sure that the code_units struct is in always in a proper state. What if we let Grackle directly track this unit information?¹

Essentially, the idea would involve 3 main changes:

• explicitly store a copy of the code_data instance used for initializing grackle_data_storage inside the instance of grackle_data_storage.
• introduce a utility function to query the units at any given time (in case the developers of the downstream application want to verify that they are evolving their unit-system consistently with Grackle's requirements)

  double * gr_required_units(chemistry_data_storage*, const char* unit_name, double current_a_value);

  • This function would return the required value of the code-units named unit_name at the specified current_a_value
  • This function would report errors (e.g. an unknown unit_name was specified or current_a_value was specified with an invalid value) by returning a negative value.
• Introduce a current_a_value member to the grackle_field_data struct. When Grackle is taking responsibility for tracking its code-units, the value stored by the member would be used to determine the current code-units during function calls to {local_}solve_chemistry or functions like {local_}calculate_temperature

I think we would need to take some care in 2 regards:

1. We would need to adopt consistent handling for current_a_value (whether it's an argument of gr_required_units or a member of the grackle_field_data struct. I would propose that a value of -1 indicates that the a_value value used during initialization should be used. I would also propose for comoving_coordinates == 0, that we treat any current_a_value value other that does not exactly match -1 or the a_value value used during initialization as an error.
2. We'd need to decide on a backwards compatible way to enable this behavior in functions like {local_}solve_chemistry or functions like {local_}calculate_temperature (i.e. when do we honor the current_a_value value stored within the grackle_field_data struct?). I see a few options:
   • In an ideal, we would just require that when agrackle_field_data is initialized with a current_a_value that disables this behavior by default and that requires users to opt into this new behavior (unfortunately, this is not a backwards-compatible option since we don't currently require users to make use of an initializer-function)
   • I think the best solution involves activating this new functionality when users explicitly pass a NULL pointer as the code_units* arg (I'm a little concerned about cases where NULL pointers are passed unintentionally, but we don't currently handle that case gracefully anyways).
   • we could potentially add a switch to the chemistry_data struct to opt into the new behavior (I'm less thrilled with this)
   • we could potentially introduce new versions of the functions that are explicitly missing the code_units struct. But, if we do that, I think it would be worth discussing some larger changes (that might eventually drop the distinction between the global API and local API).

Footnotes

1. This exact sort of automatic unit-tracking is also being suggested as part of the a proposed update to the pygrackle api in PR Proposal: More Explicit/User-Friendly Python Bindings #195. ↩

[mabruzzo]

I really want to leverage the machinery introduced by #209 for auto-computing units sooner rather than later.
To briefly recap:

• in issue Discussing Grackle's unit invariants #216, we have established that when Grackle is configured to support a comoving-unit system, the unit-system has very precise requirements about how the values of the code_units struct can change with respect to scale-factor. In fact, the values of the code_units struct are fully parameterized by the initial choice of code_units and the scale-factor (there are NO degrees of freedom)
• PR Introduce gr_required_units #209, introduces a function called gr_required_units that computes the code-units that Grackle requires at a given scale-factor. That PR is approved -- the merger process is just being held up by a very tangential testing-issue

Currently, gr_required_units is just provided as a user-convenience (to let developers confirm that they are passing in the correct value). But, I would like to use this functionality more broadly in order to simplify the API.

Essentially I want to be able to support passing a NULL pointer in place of a code_units pointer to the API functions that perform chemistry-calculations.

• These API functions include local_solve_chemistry or solve_chemistry and to functions like calculate_temperature or local_calculate_temperature.
• When passed a NULL pointer, these functions would internally compute the current units based on the current value of the scale factor.

The challenge is how do we specify the current scale-factor? I see 2 approaches:

1. We add a field called current_a_value to the grackle_field_data type
   • there is some logical sense to doing this. The values of each fluid field are only meaningful when you know the associated units (which determined by the current redshift).
   • Plus, there are scenarios where you could imagine different sets of fields might have different redshifts (think adaptive time-stepping in an AMR simulation or higher-order time integration). But this latter point is only made for arguing that they are logically consistent
   • This makes some logistical sense
2. It could also make more sense to store this information in grackle_data_storage.
   • we would probably introduce a function that looks like the following to set it

     int gr_set_current_a_value(grackle_data_storage* my_rates, double current_a_value);

• When the above function is called, we could cache the updated units (and redshift related quantities like UV-background rates) so that we don't need to recompute things in every API functions that perform chemistry-calculations.

From an API perspective, I kinda prefer choice 1 since it avoids adding "state" to the grackle-solver. With that said, I don't think adding this "state" is really a problem for any particular applications. I definitely welcome opinions! There are currently some internal advantages to state 2, but those may go away after the ongoing refactoring efforts are complete.

In all cases, the requirements and handling of current_a_value in scenarios when Grackle is configured without a comoving unit-systems would mirror the requirements of the current_a_value argument that gets passed to gr_required_units.

[gregbryan]

I like this idea a lot. It will simplify things for users and help emphasize that our choice of comoving units is a particular choice and is not (currently) super flexible. I think a neat eventual extension that Matthew mentioned earlier is to support other comoving coordinate choices. In terms of how to specify the scale-factor, I don't have a strong preference. Part of the grackle_field_data seems reasonable to me.