teos-10 / gsw-c Goto Github PK

TEOS-10/GibbsSeaWater.jl#18 (comment) and https://discourse.julialang.org/t/a-weird-issue-with-julia-1-7-1-on-ubuntu-gibbsseawater-jl-package/73632 show the error for lon > 360. The gsw_saar function is only handling negative longitudes, not longitudes beyond 360. Likely fix: change

        if (lon  <  0.0)
            lon += 360.0;

to

        lon = fmod(lon, 360);
        if (lon  <  0.0)
            lon += 360.0;

I would be inclined to fix it first in Python and then propagate it back to here, since testing is so much easier in Python.

See TEOS-10/GSW-Python#40 (comment)

I have recently tested this function using the simple examples provided at http://www.teos-10.org/pubs/gsw/html/gsw_geo_strf_dyn_height_pc.html.

I found that the returned values are not as expected.

Looking at the code, it seems that line 60 should be changed form:
prv_dyn_height_deep = (i==0)? -delta_h[0]: dyn_height_deep;
to
prv_dyn_height_deep = (i==0)? 0: dyn_height_deep;

Currently, the first value of dyn_height_deep used in the computation is -delta_h[0]*2 (see line 61) so that the resulting geo_strf_dyn_height_pc are all offset by -delta_h[0].

What is the expected behavior? Should I replace NaN values in the input data with something else?

The following test program is inspired by TEOS-10/GSW-Python#21, in which erroneous results, and sometimes memory corruption, are appearing in the GSW-Python. The test case is very simple: uniform salinity and temperature, and uniform pressure increments. In the C test program below, supply the number of points, the minimum pressure, and the pressure increment as command-line arguments. Different combinations yield different oddities. Here is an example (on OS X):

(test) efiring@manini2:~/work/programs/TEOS-10/github/c_tests$ ./gsw_test 10 0 0.5
pressure, salinity, temp., dyn height
 0.00  35.0  10.0     0.000
 0.50  35.0  10.0    -0.006
 1.00  35.0  10.0    -0.012
 1.50  35.0  10.0    -0.018
 2.00  35.0  10.0    -0.024
 2.50  35.0  10.0    -0.030
 3.00  35.0  10.0    -0.036
 3.50  35.0  10.0    -0.043
 4.00  35.0  10.0    -0.049
 4.50  35.0  10.0    -0.049

--the end--

The anomaly here is that the last two dynamic heights are identical. Now, here is more obvious problem:

(test) efiring@manini2:~/work/programs/TEOS-10/github/c_tests$ ./gsw_test 10 1 0.1
pressure, salinity, temp., dyn height
 1.00  35.0  10.0    -0.012
 1.10  35.0  10.0    -0.013
 1.20  35.0  10.0    -0.015
 1.30  35.0  10.0    -4.630
 1.40  35.0  10.0    -4.630
 1.50  35.0  10.0    -4.630
 1.60  35.0  10.0    -4.630
 1.70  35.0  10.0    -4.224
 1.80  35.0  10.0    -4.630
 1.90  35.0  10.0    -4.630

And here is the test code:

#include <stdio.h>
#include <stdlib.h>
#include "gswteos-10.h"

int main(int argc, char **argv)
{
    double *sa, *ct, *p, *dh;
    int n;
    double p_ref=0.0;
    double min_p, del_p;
    int i;

    if (argc != 4)
    {
        printf("Usage: gsw_test np min_p del_p\n");
        printf("(number of pressures, minimum pressure, pressure increment)\n");
        exit(0);
    }
    n = atoi(argv[1]);
    min_p = atof(argv[2]);
    del_p = atof(argv[3]);

    sa = malloc(n*sizeof(double));
    ct = malloc(n*sizeof(double));
    p = malloc(n*sizeof(double));
    dh = malloc(n*sizeof(double));

    for (i=0; i<n; i++)
    {
        sa[i] = 35;
        ct[i] = 10;
        p[i] = i * del_p + min_p;
        dh[i] = 0.0/0.0;  /* NaN; ensure the initial value doesn't matter */
    }

    dh = gsw_geo_strf_dyn_height(sa, ct, p, p_ref, n, dh);

    if (dh == NULL)
    {
        printf("dh is NULL\n");
        exit(-1);
    }
    printf("pressure, salinity, temp., dyn height\n");
    for (i=0; i<n; i++)
    {
        printf("%5.2f  %4.1f  %4.1f   %7.3f\n", p[i], sa[i], ct[i], dh[i]);
    }
    printf("\n--the end--\n");

    free(sa);
    free(ct);
    free(p);
    free(dh);
    return 0;
}

From @PaulMBarker in #27

Hi,

I can not work out how to start a new topic, so I am going from here.

We have a computer programmer who is going to spend some time translating the matlab code to C.
I am guessing they will not get all of the missing codes converted but they will make a contribution. I have been asked to make a list of which programmes they want done in order of importance. Do you have a preference for which programmes are prioritised ?

Paul.

If I am reading the C code correctly, gsw_internal_energy_t_exact and gsw_isochoric_heat_cap_t_exact from GSW-MATLAB are not implemented. (I did not attempt to make a complete list.) Consequently, these functions do not exist in the Python and Julia packages that wrap the C code. These functions would have usage in an introductory physical oceanography course that I teach in the MIT/WHOI Joint Program.

Does it make sense to try to implement these functions here? Is there a previous issue that gives a blueprint for adding functions?

Alternatively, there are Python or Julia wrappers that could handle the difference between gsw_internal_energy_t_exact and gsw_internal_energy in those languages instead.

I am trying to find out why we have (since some time) a test failure in GibbsSeaWater.jl.

In this package, we are using the NetCDF file to get the reference values for testing the julia package using a similar approach than this python script.

It seems however, that the GSW-C library has been updated but not the reference values in the NetCDF file.

There is another python script which read data from the matlab gsw_data_v3_0.mat and a numpy file geo_strf_dyn_height.npy.

Would it not be preferable that the reference values were accessible in a programming-language neural format? (like NetCDF)

Working with the C code is more of a pain than it needs to be for at least 2 stylistic reasons:

1. The indentation is a mix of tabs and spaces, which causes Python people like me to shudder.

2. Having both the single-file, which is compiled into the library, and the separate files for each function, which are concatenated to generate the single file, makes code modification and additions unwieldy; compilation errors and warnings refer to the single file, but all editing must be done on the individual files.

Is there support for fixing either or both of these aspects? I would be inclined to keep the separate files and ditch the concatenation.

I am working on an automated wrapping of the C code in Python, in the course of which I find some C functions that could be updated in various ways to better match the matlab, and some functions that I suspect are left over from earlier versions; they have no matching m-files. The functions I have encountered so far are:

In C only? Obsolete? Delete?

• ct_freezing_exact
• pt0_cold_ice_poly
• pt_from_pot_enthalpy_ice_poly_dh
• t_freezing_exact

Extra (optional) arguments in matlab:

• pot_enthalpy_ice_freezing mismatch: ('sa', 'p'), ('SA', 'p', 'saturation_fraction')
• z_from_p mismatch: ('p', 'lat'), ('p', 'lat', 'geo_strf_dyn_height', 'sea_surface_geopotental')

I could make PRs for such changes to bring the C closer to the current matlab; would they be welcome? (I would probably not do anything about optional output arguments; in the one case that comes to mind, it is a nearly-useless validity flag that can be obtained separately.) The method to be used for handling optional arguments is open to discussion. Options include using separate functions for each signature, or using a single variadic signature.

My priority is in the basic functionality, especially the simple scalar functions; for some of the more complex functions, including most that operate on profiles, I think that re-implementation in Python will be adequate for my purposes.

See TEOS-10/GSW-Matlab@38c9635 in the Matlab repo. This should be updated through downstream repos.

See TEOS-10/GSW-Python#48.
Are there any strenuous objections to adding the geo_strf_dyn_height and sea_surface_geopotential arguments to the p_from_z function and its inverse? They would be required arguments, so any present code using the old signatures would have to be updated. I suspect that in practice this is just the R version. @dankelley?
Although the arguments would be required in C, I would make them optional in Python.

Does GSW-C have a DOI? I would like to be able to properly cite it. The same DOI could be used to link GSW-Python, GSW-R, and others to it.

Zendo is probably the most convenient way of doing that. I would be happy to help.

Presently, line 7813 of gsw_oceanographic_toolbox.c (link: https://github.com/TEOS-10/GSW-C/blob/master/gsw_oceanographic_toolbox.c#L7813) reads

	        if (p_i[i] == p[j]) {

but comparing floating-point numbers is problematic, owing to binary representation, and I think this would be better as (e.g. judging sameness as being within 1mm in depth)

               if (fabs(p_i[i] - p[j]) < 0.001) {

To whoever is maintaining the C (and possibly fortran and matlab versions):

I ran into this because I was getting NaN values in my GSW-R code, which uses a copy of gsw_oceanographic_toolbox.c. Making the change removes the NaN problems.

I had a very quick look at the fortran code, but things are different enough there that I could not see whether the error was repeated there, also.

I did a little grepping for similar things and, ignoring cases that are at a glance about integer comparisons,

git grep -n " == " src/gsw_oceanographic_toolbox.c

gives as follows. My head is not inside the code, so I don't know if these are problems also.

src/gsw_oceanographic_toolbox.c:3717:	    if (p[ibottle] == p_ref) {
src/gsw_oceanographic_toolbox.c:3722:	if ((dp_max <= max_dp_i) && (p[0] == 0.0) && (ipref >= 0)) {
src/gsw_oceanographic_toolbox.c:3841:	            if (p[ibottle] == p_ref) ibpr = p_cnt-1;
src/gsw_oceanographic_toolbox.c:3865:	        if (p[nz-1] == p_ref) ibpr = p_cnt-1;
src/gsw_oceanographic_toolbox.c:10395:	    if (z == x[k+1])

PS. my concern about whether this problem is reproduced in some of the other comparisons, and my hope that someone might look also at the fortran and matlab code for similar problems, explains why I'm providing these notes in an issue report, as opposed to making a fork and then a pull request.

Best regards, Dan.

Question from TEOS-10/GSW-Python#3:
@dankelley or anyone else, are you building the C library on Windows? If so, are you applying any patches to do so? Or any other tricks?

Hello.

I applied the library in Arduino, more specifically the function "gsw_ps_from_c.c". However, when I compiled it it gave the following error:
C:\Program Files (x86)\Arduino\libraries\GSW-C-master/gsw_check_data.c:10148:1: warning: floating constant exceeds range of 'double' [-Woverflow]

Can anybody help me?

GSW-C currently only uses scalar inputs. During the #29, the array inputs were discussed. Maybe, I can work on vectorization using an external library (OpenBLAS). Vectorization is mostly used for optimization reasons, OpenBLAS will be a dependency but the library will be more optimized. I think I can add it as an optional dependency, so anyone who wants to use the vectorized version should compile with OpenBLAS.

So compilation of GSW-C can create two different shared objects,

• libgswteos-10.so/.dll
• libgswteos-10-vector.so/.dll

and two different API libraries,

• gswteos-10.h
• gswteos-10-vector.h

with functions such as,

void gsw_alpha(double *sa, double *ct, double *p, double *output, const uint64_t length)

It will take quite some time so would you consider adding such an feature and external (optional) dependency to the library? Additionally I'm not familiar with other implementations of GSW (for example Fortran or MATLAB etc), is there any problem/conflict with other implementations?

Assuming TEOS-10/GSW-Python#63 works, a similar wholesale update of test data
and check values will be needed here, to keep things in sync and to prevent the failure on some architectures
described in TEOS-10/GSW-Python#40.

The current Matlab code handles the following (ns, nt, np) argument sets to the gibbs function, that are not included in GSW-C:
030, 021, 111, 201, 120, 210, 102, 012.
Evidently these are not use in any of the other functions included in GSW-C.
Is there any point in adding them?

If #39 goes in, it will be helpful to have a little bit of CI. It would be enough to build and run gsw_check on Windows and Linux.

These functions (at https://github.com/TEOS-10/GSW-C/blob/master/gsw_oceanographic_toolbox.c#L5356 and https://github.com/TEOS-10/GSW-C/blob/master/gsw_oceanographic_toolbox.c#L5401) are the same, apart from a comment line and a space in one code line. Also, the matlab webpages (at http://www.teos-10.org/pubs/gsw/html/gsw_melting_ice_equilibrium_SA_CT_ratio.html and http://www.teos-10.org/pubs/gsw/html/gsw_melting_seaice_equilibrium_SA_CT_ratio.html) appear to be the same (e.g. the check values seem to match, by eye).

I think the C is just inheriting the twinning from the matlab. I assume the Fortran has the same issue. I've not checked on python, haskell, etc.

Dan.