Revising custom reference framework

[illorenzo7]

Hi all (in particular @tukss @lkorre @jorafb @feathern; please let me know if there are others who use the custom reference framework whom I may have missed). This is to follow up on a discussion some of us had over a year ago regarding the current custom reference state framework. It's been a while, so do correct me if I get anything wrong!

Just to recap, the complete reference state in Rayleigh fully determines the way the equation sets are solved, and consists of various constants and functions, as described here (https://rayleigh-documentation.readthedocs.io/en/latest/doc/source/User_Guide/physics_math_overview.html#the-system-of-equations-solved-in-rayleigh).

To be pedantic, there are also physical controls that turn various terms in the equations on and off. Most of these (with the exception of heating_type) are in Physics/Controls.F90.

A year or so ago, we determined primarily that

(A) The core functionality of the reference state interface was fine and works (mostly) as intended. This core functionality is:

1. User specifies built-in Rayleigh reference states/diffusion profiles/heating types: Set reference_type = 1, 2, or 3 and heating_type = 0, 1, or 4, and [diffusivity]_type = 1 or 2, etc. Nothing here needs modification.
2. User specifies a custom reference state contained in a binary file that specifies all relevant constants and functions (custom_reference_file = [your_file_name]). To do this, set reference_type = 4, all [diffusivity]_type's = 3, and heating_type = [any integer greater than 0].

This core functionality is fine and needs no (or very minor, for 2) updating!

The other conclusion was

(B) When the user wants more complicated things (i.e., mix-and-match reference_type = 4 with other types), it's unclear how the user does this, and indeed exactly what happens to the reference state when these mix-and-match situations are attempted (e.g., what happens if the user specifies reference_type = 4, heating_integral = 2.0d0, override_constant(10) = .true., ra_constants(10) = 3.0d0? Who gets precedence and does the code warn the user?)

I would additionally add a

(C) In the current framework, Rayleigh doesn't always throw errors when it should, or throws errors when it shouldn't. Some of the errors also don't result in the code stopping, which possibly they should? (I don't know if in Rayleigh we use warnings vs. errors, etc.)

Our goal was to update the interface (and the documentation) to address (B) and (C) above, while breaking as little that was done in the past as possible. I was in charge of coding that, and clearly did not, being distracted by a multitude of other things (namely thesis writing)! So I propose (now that I have time to look into this) to once again try to update the interface. I am reaching out to everyone above since you are the people I know use custom reference states, and the goal is still to maintain as much backward compatibility as possible.

The proposal from a year ago was:

Make the reference_type, [diffusivity]_type's, and heating_type distinct, i.e.

nu_type, kappa_type, eta_type = 3 set: (c_5, f_3, and f_11), (c_6, f_5, and f_12), and (c_7, f_7, and f_13), respectively (and treat the new scalar diffusion coefficients similarly)
heating_type = 5 (a new type) sets: c_10 and f_6
reference_type = 4 sets: everything else

Then the mix-and-match choices (e.g., reference_type = 4, heating_type = 1) would yield (e.g.) most of the functions/constants from the custom_reference file, but the heating Q(r) proportional to rhoT = f_1f_4.

Note that mix-and-match between reference_type = 4 and the [diffusion]_type's DOES currently work as described above and I don't think needs changing. Thus, I believe the only functionality that people used that we agreed to change would be to add a new

heating_type = 5

The only use case that would change that I can see would be, for example

reference_type = 4
heating_type = 1

Before, that would read c_10 and f_6 from the reference file. Under the new proposal, heating would be set based on heating_type = 1: i.e., proportional to rho*T.

There is also the (I think separate) question of the various modifications to the constants. These come from override_constants, with_custom_reference, etc. I think the only thing that needs to be done there is to specify a hierarchy. I'll make a separate issue for that.

For now, I propose only:

add heating_type = 5 (to read the heating from the custom reference file)

and furthermore, allow heating_type's = 1, 4 to be respected when reference_type = 4.

From what I remember, this would break a few cases only for @lkorre and me (@illorenzo7), and we were willing to make the sacrifice. Please let me know if you use the custom reference stuff and foresee problems! Otherwise, I'd like to get to work on a pull request and have something for people to review fairly soon.

-Loren