Re-factor some EOS functions to not directly port GSW-C Toolbox code (license violation)

[katsmith133]

The GSW-C license permits redistribution "in source and binary forms, without modification". Porting the code into a kokkos kernel is modification (thank you to @xylar for pointing this out). So all TEOS-10 functions we implement need to be derived in some form from Roquet et al. 2015, the original 75-term polynomial for the EOS, and not referenced directly from the GSW-C toolbox code. The currently implemented functions calcGibbsDerivPt0Pt0, calcPtFromCt, calcCtFromPt, and calcCtFreezing are all likely in violation of this license's terms, and need to be re-factored similar to the specific volume derivatives @xylar has implemented in #485.

[xylar]

@katsmith133 and @alicebarthel, I have an idea on how to get around the gsw-c license.

First off, Claude tells me that the full TEOS-10 implementations are all published. So those would always be available as a fallback.

But we could also do our own polynomial approximation of any functions that aren't published. That's the sticking point. Some TEOS-10 functions have polynomial approximations that are not published outside of gsw-c. But we could make our own.

[xylar]

Here's an early draft of that idea:
https://github.com/xylar/E3SM/blob/omega/design/teos-10-emulator/components/omega/doc/design/Teos10Fits.md

[alicebarthel]

@xylar thanks for thinking about this. The early teos-10 was documented in the paper so I had no issue with it but I am not sure what Claude means saying the rest has been published (beyond the gsw-c code base). I have reached out to get guidance on options to contribute code back to the repo. I'll have a look at your suggestion but happy to meet online to discuss.

[xylar]

@alicebarthel, happy to chat more the next time we meet.

Response from Claude:

That's a fair challenge, and my wording above was too loose. Let me be specific, because there are two different things that could be meant by "published" and only one of them is fully public.

1. The standard itself — published, and enough on its own.

The TEOS-10 primitives are freely available releases, not code:

• Seawater Gibbs function: IAPWS R13-08 (PDF mirrored at teos-10.org)
• Ice Ih Gibbs function: IAPWS R10-06(2009), and Feistel & Wagner (2006), J. Phys. Chem. Ref. Data 35, 1021–1047
• TEOS-10 Manual (IOC/SCOR/IAPSO 2010, Manuals and Guides 56)
• The frazil/melting thermodynamics: McDougall, Barker, Feistel & Galton-Fenzi (2014), JPO 44, 1751–1775

Every function we need is defined by these, with full coefficient tables, and none of them require reading GSW-C. That's the part I meant, and I should have said so.

2. GSW's polynomial approximations to the standard — only partly published.

This is the part you're right to doubt, and where I overstated. Roquet et al. (2015) published theirs, which is why our EOS is fine. Most of the others exist only as coefficients in the code.

The frazil function is a good worked example. TEOS-10 publishes notes on gsw_frazil_properties_potential_poly, which derive the whole method — the bulk enthalpy relation, the func0 test, the Newton scheme — and describe the approximation exactly:

In the GSW code we first estimate the ice mass fraction using a polynomial with 36 coefficients. This polynomial is in terms of (normalized values of) SA_B, func0 and p. Every term in the polynomial contains at least the first power of func0 so that w_Ih → 0 as func0 → 0.

So the form, the variables, the continuity constraint, and even the fact that they upweighted the parameter range ocean models actually use are all published. The 36 numbers are not. Same for the 29-coefficient derivative polynomial described in the same document, and gsw_pot_enthalpy_ice_freezing_poly cites only the TEOS-10 Manual and a Newton-method paper for its coefficients.

Why this doesn't block us: the unpublished layer is exactly the layer we don't need. We can't copy those coefficients, but we don't have to — layer 1 defines the answer, so we can fit our own approximations to it and use the unmodified toolbox as an oracle to check them, which the license explicitly permits. The publication status of GSW's polynomials only sets how much work we do, never whether we're stuck.

On the existing functions. @katsmith133's original point checks out concretely, at least for calcCtFreezing: its coefficients C0–C22, A and B are the values of GSW_FREEZING_POLY_COEFFICIENTS in GSW-C's gsw_internal_const.h, transcribed verbatim under renamed symbols. I've updated the design doc accordingly — I had originally (and wrongly) treated calcCtFreezing as a safe building block for the frazil work, when in fact it's on the critical path and has to be re-derived first.

On relicensing. Very glad you're pursuing that, and it's the cleanest outcome if it lands. Worth noting it would shrink this work rather than remove it: it would let us port the cheap polynomial functions directly, but the expensive ones are expensive because of iteration, not arithmetic. gsw_pt_from_t_ice runs up to six Newton steps on the ice Gibbs function and the frazil solve runs three more, each needing the freezing temperature and ice freezing enthalpy. Per-thread-variable trip counts are bad on GPU regardless of who owns the copyright, so we'd still want fixed-cost replacements for those. So the two paths are complementary, and neither is blocked on the other.

The design doc is here; sections 3.1 and 3.1.1 are the relevant ones, and §6 now has the full reference list.

[xylar]

Written by Claude Code and posted from @xylar's account. The Aurora profiling summarized here was also run by an agent; happy to share the full results and caveats.

Update: we profiled the performance premise before building anything, and it didn't hold, so we've changed approach.

The earlier draft proposed a general polynomial-fitting framework, motivated partly by licensing and partly by an expectation that fitted approximations would be cheaper than GSW's algorithms on GPU. Measured on Aurora (one PVC tile, fp64, Omega's own launch geometry):

• Teos10Eos is memory-bandwidth-bound — 1.06–1.26× a pure streaming kernel with the same memory traffic, once the working set leaves L2. Its arithmetic is essentially hidden, so a cheaper polynomial buys nothing.
• The frazil equilibrium solve, which the framework was built around, gains only 5–7% from a fit in a realistic sparse ice field. (Measured in the EOS kernel's launch geometry; add a teos-10 version of frazil formation+melt #462 maps lanes to whole columns instead, so its divergence pattern differs and this is worth re-checking there.) Divergence does serialize on the hardware, but GSW varies its trip count only between 1 and 3 and the ordinary ocean case is 1, so there's little to act on.
• One exception: gsw_pot_enthalpy_from_pt_ice_poly measured 4.4×, and not for divergence reasons — GSW Newton-inverts its own degree-7 polynomial through 5 unconditional iterations, ~16 polynomial evaluations for a 1-D function. Fitting the inverse directly is shorter and faster.

One worthwhile fit doesn't justify a fitting framework, so it's gone. The new design is licence-driven: implement each function directly from the published standards (IAPWS R13-08, IAPWS R10-06, the TEOS-10 Manual, McDougall et al. 2014), verify against the unmodified toolbox as an oracle, and fit only where a fit is simpler or measurably faster in the kernel it actually runs in.

Draft: Teos10Functions.md

Two practical notes. calcCtFreezing is first in the queue — it's one of the four functions here, it's on the critical path for the frazil equilibrium test, and #462 already calls it for every level of every cell. And none of this depends on the relicensing outcome, @alicebarthel; if that lands it shrinks the work further rather than changing its direction.