marbl-ecosys / MARBL

Marine Biogeochemistry Library
https://marbl-ecosys.github.io
Other
14 stars 25 forks source link

Refactor computation of K #470

Closed mnlevy1981 closed 1 month ago

mnlevy1981 commented 2 months ago

Paulot et al define

K = (1/kg + H/kw)^-1

rewriting this expression as

K = ((kw + kg*H) / (kg*kw))^-1 = (kg*kw) / (kw + kg*H)

allows for situations where either kg or kw (but not both) is zero. It's also computationally more efficient, replacing three divisions and an addition with two multiplications, an addition, and a division (effectively replacing two divisions with multiplications)

mnlevy1981 commented 2 months ago

Testing:

The stand-alone test suite reports round-off level changes in a small number of variables:

$ ./netcdf_comparison.py -b ../tests/regression_tests/call_compute_subroutines/baseline.nc -n ../tests/regression_tests/call_compute_subroutines/history_1inst.nc --strict exact
Comparing ../tests/regression_tests/call_compute_subroutines/history_1inst.nc to the baseline ../tests/regression_tests/call_compute_subroutines/baseline.nc

Variable: NHx_SURFACE_EMIS
... Values are not the same everywhere
    Max relative error: 2.1116423442199165e-16
    Max absolute error: 4.235164736271502e-22

Variable: STF_ALK
... Values are not the same everywhere
    Max relative error: 5.319968357761115e-16
    Max absolute error: 4.235164736271502e-22

Variable: STF_ALK_ALT_CO2
... Values are not the same everywhere
    Max relative error: 5.319968357761115e-16
    Max absolute error: 4.235164736271502e-22

Variable: STF_NH4
... Values are not the same everywhere
    Max relative error: 2.862036805523897e-16
    Max absolute error: 4.235164736271502e-22

Differences found between files!

I'll run a test through the CESM test suite as well, but the stand-alone test results matched expectations (this is a small round-off level change in an intermediate computation)

mnlevy1981 commented 2 months ago

I created a CESM test baseline with SMS_Ld1.TL319_t232.C1850MARBL_JRA, and then compared the proposed fix here and got the [expected] small differences; run $ grep RMS /glade/derecho/scratch/mlevy/tests/one-off/marbl_updates/SMS_Ld1.TL319_t232.C1850MARBL_JRA.derecho_intel.C.marbl_update/SMS_Ld1.TL319_t232.C1850MARBL_JRA.derecho_intel.C.marbl_update.mom6.h.bgc.z.0001-01.nc.cprnc.out for details, but the biggest differences reported by cprnc is

$ grep RMS /glade/derecho/scratch/mlevy/tests/one-off/marbl_updates/SMS_Ld1.TL319_t232.C1850MARBL_JRA.derecho_intel.C.marbl_update/SMS_Ld1.TL319_t232.C1850MARBL_JRA.derecho_intel.C.marbl_update.mom6.h.bgc.z.0001-01.nc.cprnc.out | grep -v E-1[6789]$ | grep -v E-2[01]$
 RMS Fe_scavenge                      3.0009E-27            NORMALIZED  1.2915E-15
 RMS Lig_scavenge                     3.4627E-29            NORMALIZED  2.2027E-15

I think this confirms our expectation; this change is not bit-for-bit, but it's definitely round-off level.