Shielding the user/coder from the "ref" object

[illorenzo7]

I am tagging @feathern @cianwilson, and @BWHindman on this because I feel they might have some ideas. A lot of the issues that have come up recently related to documentation of the diagnostics (e.g. PR #626 and #627), places where the diagnostics are physically inconsistent across reference_types (e.g. issue #628), difficulty adding new diagnostics (e.g. PR #625), and the difficulty now modifying PDE_Coefficients.F90 (just my opinion!) comes from the following:

The ReferenceInfo-type "ref" object (and nu/kappa/etc) in Rayleigh is internal to the code, but the modern user interface is purely in terms of ra_constants and ra_functions. In order to check diagnostics and make new reference_types or diagnostics, the user/coder has to constantly interpolate between ra_constants and ra_functions and the ReferenceInfo object "ref"/nu/kappa/etc. This is time-consuming, error-prone, and, in my opinion, has clearly causing been causing some real issues with adding new diagnostics, modifying PDE_Coefficients.F90, and making sure the documented diagnostics are actually what are being computed.

I have some proposals to fix this problem (which, to be clear, is not a bug, just something that seems to unnecessarily make our lives more difficult!) These would be a series of small pull requests which at NO STAGE would modify the functionality of the code at all, except at the very end for fixing some physically inconsistent diagnostics. These pull requests would slowly shield the user from having to translate between the ref object and ra_functions/constants again.

1. First, in PDE_Coefficients.F90, make a completely redundant function Set_Complete_Ref() which would set ref, nu, kappa, etc. from the available ra_constants and ra_functions. This would be called after Set_Reference_Equation_Coefficients and Set_Diffusivity_Equation_Coefficients (which do the reverse, setting ra_constants and ra_functions from the ref/nu/kappa objects), and so should do nothing at first. This function effectively would effectively be an expanded version of Get_Custom_Reference(), but set nu/kappa/etc as well.
   2-8: One by one, for subroutines corresponding to Reference_Types = 1, 2, 3, 5 and then the diffusivities and heating (respectively: Constant_Reference(), Polytropic_Reference(), Polytropic_ReferenceND(), Polytropic_ReferenceND_General(), and Initialize_Diffusivity(), and Initialize_Reference_Heating()), make the following change: initialize the appropriate ra_constants and ra_functions, in parallel with the ref/nu/kappa objects. Again this would do nothing but add redundancy. Note that reference_type = 4 doesn't need to be touched here because it works exactly by setting ra_constants and ra_functions to begin with. Note that this would effectively undo some work I did a few years ago to remove all these ra_constants/ra_functions settings from the various subroutines :)
2. Remove the now completely redundant functions Set_Reference_Equation_Coefficients() and Set_Diffusivity_Equation_Coefficients(), and calls to them. (I.e., completely undo, but not break, my previous work in PR Safeguard equation_coefficients file to reflect the actual code #447)
   10-16: One by one, remove all places the ref/nu/kappa objects are set in the Reference_Types=1,2,3,5/diffusivity/heating subroutines. These would now be automatically set by the new routine Set_Complete_Ref, and so wouldn't need to be set by the subroutines.
3. For the diagnostics (but perhaps for all modules?) we make the arrays ra_functions/ra_constants arrays visible, in addition to the ref/nu/kappa objects.
4. In various diagnostics (like the ones mentioned in issue Making diagnostic codes work across reference_types? #628) we begin replacing calls to ref%density, ref%dpdr_w_term, pressure_specific_heat, etc. by the much more readable (and always physically consistent) ra_constants(3) and ra_constants(2), etc.
5. @cianwilson I would have to chat with you about how to amalgamate the new scalar fields stuff into all this.

I understand that this would be a large undertaking, but I am happy to do it over the next few months, if others might be willing to review it. I think this approach has several advantages:
a. Nothing about how the code runs actually changes (for any of the many pull requests), except, in 18, for making certain diagnostics physically consistent across reference_types.
b. If someone wants to make new diagnostics, it will be much easier to just write down what they want in the documentation (like c_2f_2tvar*u_r) and then write basically the exact same thing in the Fortran diagnostics files. This will ensure that future diagnostics are both physically consistent across reference_types and consistent with the documentation.
c. If the user wants to create a new reference_type (like a new nondimensionalization), it will be far easier, because all they have to do is write down the equation set they want to solve in the documentation in terms of c's and f's, and then set the ra_constants and ra_functions directly in the new reference_type. This will ensure physical consistency between the documentation and the new reference type automatically, and the user won't ever have to worry about the ref/nu/kappa objects because they are automatically taken care of by the universal function Set_Complete_Ref().
d. If future bugs/inconsistencies ARE found in the diagnostics, they will be much easier to track! The subroutine Set_Complete_Ref() will also act as an easily-readable "dictionary" so the user CAN translate between ra_constants/ra_functions and "ref" when necessary and identify bugs in the diagnostics. Then they can improve code readability (and fix bugs at the same time) by replacing the calls to ref/nu/kappa with the appropriate ra_constants/ra_functions in the diagnostics Fortran code.

I know it's a lot, but thanks for reading! I really think this could make things easier for us in the future. If someone will review these many pull requests (thinking of @cianwilson specifically!), I'll do the core work :)

[cianwilson]

I would support this effort (and, in fact, would suggest we go further than the diagnostics and apply it to the equations too - though this might come at the expense of some extra flops)! In trying to work out the curl of the momentum diagnostics it was really annoying to have to translate from the custom reference state documentation to what was actually in the code. It would have been much easier to do as @illorenzo7 suggests and just write down and expand upon what was in the docs.