Physical constants

[pbosler]

Question

What would be the best (or at least acceptable to a majority) way to handle physical constants in a consistent manner throughout scream?

This issue has come up in several recent discussions (e.g., @PeterCaldwell, @AaronDonahue ) and at least 2 issues:
haero issue #6 ( @jeff-cohere and @pbosler )
scream pr #579 (@tcclevenger, @ambrad)

For discussion:

1. Currently, different parameterizations often define the same constants, but with different numbers of significant figures and/or units. The ideal solution is to use NIST-defined values for each constant, defined in a single location (file or struct) that dynamics and parameterization can use/alias with local variables as needed.
2. Suppose people agree with that statement; where should such code live?
   Physical constants have nothing to do with Kokkos, but EKAT already contains other supporting features like ekat::units that are also unrelated to performance portability, but generally useful to everything in E3SM, so I'm inclined to put physical constants in EKAT (@bartgol).
3. We've found that some constants require double precision (e.g., they may be very small or very large, or they may be required by possibly ill-conditioned algorithms) whereas others can be single precision. How should these requirements be handled in scream & ekat?
4. What should the timeline be for this work?
5. Should this be a forward looking c++ only change (to avoid BFB issues), or should it include fortran as well?

[ambrad]

I'll take the easy one, 3. Just define it as hard-coded double rather than Real. Do that consistently in F90 and C++ when needing BFB tests. If implicit type conversion in the presence of multiple types is not the same in C++ and F90 by default, I believe you can add an F90 flag or flags to get it to be the same.

[pbosler]

Oof. BFB. Added item 5 to the discussion.

[PeterCaldwell]

I agree this stuff should live in EKAT. I agree with Andrew about making params double by default. I don't think we can be "forward thinking" by just doing this in C++ because of our BFB tests. I'd suggest just copying whatever E3SM uses for our initial version of this constants database (rather than matching NIST) in order to circumvent the BFB issue while still targeting C++ only. Then when we jettison F90, we can update our constants file with NIST values.

One problem with what I just said is that E3SM may use different values for the same constant in different places in the code. And collecting E3SM values may be a hassle.

[pbosler]

E3SM may use different values for the same constant in different places in the code.

It does. This means that there will be BFB problems no matter what we do.

And collecting E3SM values may be a hassle.

It is, but I've already done this for MAM. It's not a particularly exciting day at work, but it's not too bad with some decent background music.

[bartgol]

Yes, Ekat already contains stuff that has nothing to do with Kokkos: yaml, parameter lists, units, std meta-utils,...

I actually pondered on whether units should be in Ekat. Eventually I left them there, but I don't feel strongly either way.

Regarding constants, I guess Ekat can be a reasonable place to put them. If they have to be used in the atm only, then scream/src/share is also a good place. Recall that Ekat has no knowledge of single precision vs double precision build, so if you store constants in ekat, you either store them all as double, or you need to put them inside a templated class:

namespace Constants1 {
static constexpr double g = 9.80665;
}

template<typename RealType>
struct Constants2 {
   static_assert(std::is_same<RealType,double>::value || std::is_same<RealType,float>::value,
                         "Error! Only 'float' or 'double' supported.\n");
   static constexpr RealType g = 9.80665;
};

Imho, the latter is superior (app-side code can always do using Constants = ekat::Constants2<double>;).

[pbosler]

put them inside a templated class

This is similar to our current approach in haero, but there we don't encapsulate the constants in a struct -- they're namespace variables like the ekat units instantiations.

static constexpr double dbl_required_const = val;
static constexpr Real flexible_type_const = other_val;

If we decide to use structs instead, this means we would have multiple structs --- one for required double precision, one that has a template parameter to allow for single precision (which is ok, it's just something the ekat clients would have to know).

edit: It's possible there are additional meta programming tricks that I'm currently unaware of.

[bartgol]

Are the constants required in DP always going to be needed in DP only? If so, then yes, you need 2 structs (one for DP only, one templated), otherwise you can just stuff everything in the templated one.

[pbosler]

I'm not sure how far this kind analysis has gone in the broader community... But here's an example:

program avogadro_test

real(4) :: sp_avo = 6.2214E26
real(8) :: dp_avo = 6.2214D26

write(6,*) "abs. diff.:           dp(N_A) - sp(N_A) = ", dp_avo - real(sp_avo,8)
write(6,*) "rel. diff.: (dp(N_A) - sp(N_A))/dp(N_A) = ", (dp_avo - real(sp_avo,8))/dp_avo

end program

Output:

$ gfortran avogad.f90 -o avo
$ ./avo
abs. diff.: dp(N_A) - sp(N_A) = 1.8068707374507491E+019
rel. diff.: (dp(N_A) - sp(N_A))/dp(N_A) = 2.9042831797517423E-008

So, a single-precision Avogadro constant represents ~2E19 fewer molecules per kmol than a double-precision version. Over every column, parameterzation, species, and model level, that adds up to a lot of molecules.

[pbosler]

@PeterCaldwell , @bartgol

Regarding constants, I guess Ekat can be a reasonable place to put them. If they have to be used in the atm only, then scream/src/share is also a good place.

This is an important point. In mpas/src/framework/mpas_constants.F, the radius of the sphere is defined to be 6371229.0 meters, but in e3sm/cime/src/share/util/shr_const_mod.F90, it's 6.37122e6. Is E3SM modeling a planet with 9 meters of sea level rise at t=0?

[ambrad]

@pbosler only relative error matters. Your point about condition number is the important one. An absolute "big" or "small" number doesn't mean much.

[pbosler]

@ambrad I understand, and the absolute difference I mention primarily to point out the inconsistencies that exist in the code already. At the same time, I'm not sure we (as in the climate modeling community "we") understand well what is or is not ill-conditioned in the coupled model.

[bartgol]

@PeterCaldwell , @bartgol

Regarding constants, I guess Ekat can be a reasonable place to put them. If they have to be used in the atm only, then scream/src/share is also a good place.

This is an important point. In mpas/src/framework/mpas_constants.F, the radius of the sphere is defined to be 6371229.0 meters, but in e3sm/cime/src/share/util/shr_const_mod.F90, it's 6.37122e6. Is E3SM modeling a planet with 9 meters of sea level rise at t=0?

I don't think anyone thinks we should pick up the challenge of making constants uniform across e3sm. Imho, scream could simply stuff them in scream/src/share for now. However, this does not help EAGLES/Haero developers. So if those constants are needed in EAGLES/Haero (and I assume they might be), we can put them in Ekat.

I feel like Ekat is becoming more and more of a E3SM std library, with features across a wide range of topics. Maybe the name should be re-considered at this point.

[pbosler]

@bartgol I agree that we can't change E3SM without significant help from the core/infrastructure teams; EKAT and SCREAM are good opportunities to demonstrate these design principles so that, perhaps, they might be adopted later. It think it's enough (and certainly challenging enough) to try for a uniform treatment throughout scream.