search

TEOS-10 / GSW-R

R implementation of the Thermodynamic Equation Of Seawater - 2010 (TEOS-10)
http://teos-10.github.io/GSW-R/
Other
8 stars 5 forks source link

gsw-3.05 should be used #14

Closed dankelley closed 8 years ago

dankelley commented 9 years ago

Goal. The C and Fortran codes for GSW-3.05 are now available on the TEOS-10 website, so it makes sense to start migration of GSW-R to 3.05. I made a new branch named 305 for this. Below is a list of the functions in gsw (from a grep in the file R/gsw.R), with a checkmark indicating that the function and its documentation have been updated, and that the check values match the published ones (but see the notes on some functions).

dankelley commented 9 years ago

I wrote some code to get the data from the netcdf file in the fortran 3.05 version. It is at https://github.com/TEOS-10/GSW-R/blob/305/create_data/create_saar.R in case @richardsc wants to check my work. The numbers in the file are the same as they were for version 3.02, at least to the 1e-7 tests that are done by default with the testthat package in R.

PaulMBarker commented 9 years ago

FYI, the netcdf file is also available on the github under datafile.

dankelley commented 9 years ago

@fdelahoyde and @PaulMBarker ...

As an update, I've updated gsw-r for the new C code. Part of that work is to use the TEOS-10 supplied check values in build-time tests. By default, in R build-tests use a relative error of sqrt() machine precision, or 1.5e-8 on my computer. The functions all check OK at this threshold, with one exception: gsw_latentheat_melting(), which has deviations about 20 times the threshold. (Details below, in funky R notation but I think it will be clear enough.)

Paul -- is there a chance the webpage http://www.teos-10.org/pubs/gsw/html/gsw_latentheat_melting.html might be wrong? (I don't have matlab so I can't check.)

Frank -- could you see if you get the same values from C as I do? (Possibly I made an editing error. I had to do a few changes in the code so that R can save the data file separately, because the R core team won't permit a C file that contains so much data.)

> library(gsw)
Loading required package: testthat
> SA <- c(34.7118, 34.8915, 35.0256, 34.8472, 34.7366, 34.7324)
> p <- c(      10,      50,     125,     250,     600,    1000)
> lh <- gsw_latentheat_melting(SA, p)
> lh_check <- 1e5*c(3.299495187300804, 3.298611954422526, 3.297124383647952,
+                   3.294972884747496, 3.288480015369891, 3.280715953443947)
> data.frame(lh=lh, difference=lh-lh_check)
        lh   difference
1 329949.7  0.149297041
2 329861.3  0.139797546
3 329712.6  0.123918659
4 329497.4  0.101058326
5 328848.0  0.043018985
6 328071.6 -0.009102756
> mean(abs((lh-lh_check)/lh))
[1] 2.862329e-07
PaulMBarker commented 9 years ago

Hey Dan,

You are correct the values listed on the webpage for latent heat are wrong. I must have missed updating the values when I updated the toolbox, sorry about that.

Here are the correct values:

SA = [34.7118; 34.8915; 35.0256; 34.8472; 34.7366; 34.7324;] p = [ 10; 50; 125; 250; 600; 1000;]

latentheat_melting = gsw_latentheat_melting(SA,p)

latentheat_melting =

1.0e+05 *

3.299496680271213 3.298613352397986 3.297125622834541 3.294973895330757 3.288480445559747 3.280715862416388

By the way when I changed the datafile to a netcdf file I included lots of test values and the computational accuracy that a calculation should achieve, these values will always be more upto date than the web page examples as they are automatically generated when I make the datafile, In the last couple of versions I have written more code to automate these things. Initially I was doing much of it by hand. Unfortunately the webpages are still done one at a time by hand.

It is a pain not getting recognition for your hard work translating GSW, I can assure you we appreciate it.

dankelley commented 9 years ago

Thanks, Paul -- the values work with the default equality-check tolerances.

dankelley commented 9 years ago

@PaulMBarker -- I wonder if the web check values for gsw_t_freezing() may be in error. I get as follows (mean absolute error a bit under 1e-4 degC) when I hook R to gsw_t_freezing(), but things are OK if I hook to gsw_t_freezing_exact().

> library(gsw)
> library(testthat)
> SA <- c(34.7118, 34.8915, 35.0256, 34.8472, 34.7366, 34.7324)
> p <- c(      10,      50,     125,     250,     600,    1000)
> saturation_fraction <- 1
> tf <- gsw_t_freezing(SA, p, saturation_fraction)
> expect_equal(tf, c(-1.902730710149803, -1.942908619287183, -2.006861069199743,
+                    -2.090985086875259, -2.351293130342102, -2.660498762776720))
Error: tf not equal to c(-1.9027307101498, -1.94290861928718, -2.00686106919974, -2.09098508687526, -2.3512931303421, -2.66049876277672)
6/6 mismatches (average diff: 7.24e-05).
First 6:
 pos     x     y      diff
   1 -1.90 -1.90 -2.63e-05
   2 -1.94 -1.94 -1.07e-05
   3 -2.01 -2.01  1.63e-05
   4 -2.09 -2.09  5.63e-05
   5 -2.35 -2.35  1.38e-04
   6 -2.66 -2.66  1.86e-04
PaulMBarker commented 9 years ago

I need a bit more info to help on this one.

SA = [34.7118, 34.8915, 35.0256, 34.8472, 34.7366, 34.7324]; p = [ 10, 50, 125, 250, 600, 1000]; saturation_fraction = 1;

t_freezing = gsw_t_freezing(SA, p, saturation_fraction)'

t_freezing =

-1.902730710149803 -1.942908619287183 -2.006861069199743 -2.090985086875259 -2.351293130342102 -2.660498762776720

t_freezing_poly = gsw_t_freezing_poly(SA, p, saturation_fraction)'

t_freezing_poly =

-1.902704434299200 -1.942897945475226 -2.006877364578649 -2.091041391538033 -2.351431560178946 -2.660685237633709

I do not have a gsw_t_freezing_exact, what does your t_freezing_exact code do, does it call the chemical potential of water and the gibbs function or does it use a polynomial ?

dankelley commented 9 years ago

Thanks, Paul. The R gets gsw_t_freezing_exact() from line 10046 of gsw_oceanographic_toolbox.c, which is linked with:

https://github.com/TEOS-10/GSW-C/blob/master/gsw_oceanographic_toolbox.c#L10046

and I think is derived from the fortran source linked with:

https://github.com/TEOS-10/GSW-Fortran/blob/master/toolbox/gsw_t_freezing_exact.f90

In answer to the question, it seems to be using chemical potential.

I'm not sure why the fortran and C have a function that's not in the matlab.

I'll add an @fdelahoyde to this comment so that Frank gets an email. Possibly this issue is relevant to him, too.

PS. when replying over email, it's best not to include the material that was sent. Actually, the best is just to click the github link at the end, which shows all the discussion, possibly with images and so forth. I often chop out content from others' comments, if it contains phone numbers etc.

PaulMBarker commented 9 years ago

ok that helps, gsw_t_freezing_exact is gsw_t_freezing, I will email Glenn. Yesterday I calculated the values that I sent you not just copying them from the webpages, so they should be correct.

I am not sure what i included in my reply, I will check in future not to include the received email. This is the first time I have used the web interface, but I really need the spell checker option.

dankelley commented 8 years ago

I realize this discussion concluded months ago, so I'm closing the issue.