From f8dbfecaa537ab0f55e626ae973cf70497c6fd8c Mon Sep 17 00:00:00 2001 From: pvthinker Date: Thu, 6 Apr 2023 16:56:29 +0200 Subject: [PATCH 1/4] Add Gibbs potential to the package --- gsw/__init__.py | 1 + gsw/gibbs_function.py | 58 +++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 59 insertions(+) create mode 100644 gsw/gibbs_function.py diff --git a/gsw/__init__.py b/gsw/__init__.py index 7d28f21..3f8d68d 100644 --- a/gsw/__init__.py +++ b/gsw/__init__.py @@ -36,6 +36,7 @@ from .stability import * from .geostrophy import * from .utility import * +from .gibbs_function import gibbs from . import geostrophy from . import utility from . import stability diff --git a/gsw/gibbs_function.py b/gsw/gibbs_function.py new file mode 100644 index 0000000..ba44e8d --- /dev/null +++ b/gsw/gibbs_function.py @@ -0,0 +1,58 @@ +""" +Wrapper to provide the Gibbs potential +""" +import numpy + +from . import _gsw_ufuncs +from ctypes import CDLL, c_double + + +def _set_gibbs(): + name = "gsw_gibbs" + libname = _gsw_ufuncs.__file__ + libobj = CDLL(libname) + gibbs_compiled = getattr(libobj, name) + gibbs_compiled.restype = c_double + + def gibbs(ns, nt, np, SA, CT, p): + """ + Calculates the specific Gibbs free energy and derivatives up to order 2 + + Parameters + ---------- + ns : int + order of SA derivative + nt : int + order of CT derivative + np : int + order of p derivative + SA : array-like + Absolute Salinity, g/kg + CT : array-like + Conservative Temperature (ITS-90), degrees C + p : array-like + Sea pressure (absolute pressure minus 10.1325 dbar), dbar + + Returns + ------- + gibbs : specific Gibbs energy [J/kg] or its derivative + + + """ + if isinstance(SA, numpy.ndarray): + shape = SA.shape + assert CT.shape == SA.shape + assert p.shape == SA.shape + value = numpy.zeros(shape) + for k in numpy.ndindex(shape): + value[k] = gibbs_compiled(ns, nt, np, + c_double(SA[k]), c_double(CT[k]), c_double(p[k])) + else: + value = gibbs_compiled(ns, nt, np, + c_double(SA), c_double(CT), c_double(p)) + + return value + return gibbs + + +gibbs = _set_gibbs() From de5a2226ae1c0adf904d927dc3831cd5e7508fd1 Mon Sep 17 00:00:00 2001 From: pvthinker Date: Thu, 6 Apr 2023 17:13:13 +0200 Subject: [PATCH 2/4] Add test to gibbs --- gsw/tests/test_gibbs.py | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) create mode 100644 gsw/tests/test_gibbs.py diff --git a/gsw/tests/test_gibbs.py b/gsw/tests/test_gibbs.py new file mode 100644 index 0000000..b43de6c --- /dev/null +++ b/gsw/tests/test_gibbs.py @@ -0,0 +1,28 @@ +import numpy as np +from numpy.testing import assert_almost_equal +from gsw import gibbs + +saref = 35. +ctref = 15. +pref = 0. + + +def test_gibbs(): + value = -1624.8303316 + assert_almost_equal(gibbs(0, 0, 0, saref, ctref, pref), value) + + +def test_array(): + value = -1624.8303316 + shape = (4, 5, 3) + sa = saref*np.ones(shape) + ct = ctref*np.ones(shape) + p = pref*np.ones(shape) + + res = gibbs(1, 1, 0, sa, ct, p) + assert res.shape == shape + + +def test_gibbsderivative(): + value = 0.590312110989 + assert_almost_equal(gibbs(1, 1, 0, saref, ctref, pref), value) From 23954a980589eae2cae3df829a139fcc750018fc Mon Sep 17 00:00:00 2001 From: pvthinker Date: Fri, 7 Apr 2023 11:48:14 +0200 Subject: [PATCH 3/4] Add ufunc-ificated gibbs function --- gsw/__init__.py | 3 ++- gsw/gibbs_function.py | 17 ++++++++++++++++- 2 files changed, 18 insertions(+), 2 deletions(-) diff --git a/gsw/__init__.py b/gsw/__init__.py index 3f8d68d..71696f7 100644 --- a/gsw/__init__.py +++ b/gsw/__init__.py @@ -36,7 +36,7 @@ from .stability import * from .geostrophy import * from .utility import * -from .gibbs_function import gibbs +from .gibbs_function import gibbs, gibbs_ufunc from . import geostrophy from . import utility from . import stability @@ -52,3 +52,4 @@ from ._version import __version__ except ImportError: __version__ = "unknown" + diff --git a/gsw/gibbs_function.py b/gsw/gibbs_function.py index ba44e8d..ca35a74 100644 --- a/gsw/gibbs_function.py +++ b/gsw/gibbs_function.py @@ -7,12 +7,20 @@ from ctypes import CDLL, c_double -def _set_gibbs(): +def _get_compiled_gibbs(): name = "gsw_gibbs" libname = _gsw_ufuncs.__file__ libobj = CDLL(libname) gibbs_compiled = getattr(libobj, name) gibbs_compiled.restype = c_double + return gibbs_compiled + + +_gibbs_compiled = _get_compiled_gibbs() + + +def _set_gibbs(): + gibbs_compiled = _get_compiled_gibbs() def gibbs(ns, nt, np, SA, CT, p): """ @@ -55,4 +63,11 @@ def gibbs(ns, nt, np, SA, CT, p): return gibbs +def _basic_gibbs(ns, nt, np, SA, CT, p): + """ this doc is lost during the ufunc-ification""" + return _gibbs_compiled(ns, nt, np, c_double(SA), c_double(CT), c_double(p)) + + gibbs = _set_gibbs() + +gibbs_ufunc = numpy.frompyfunc(_basic_gibbs, 6, 1) From 79785b733f68dafbd6158fdd75a9ee76daf41a16 Mon Sep 17 00:00:00 2001 From: pvthinker Date: Sat, 8 Apr 2023 11:26:08 +0200 Subject: [PATCH 4/4] Make `gibbs` behaves like a ufunc while having the correct doc --- gsw/__init__.py | 2 +- gsw/gibbs_function.py | 83 +++++++++++++++++++------------------------ 2 files changed, 38 insertions(+), 47 deletions(-) diff --git a/gsw/__init__.py b/gsw/__init__.py index 71696f7..e69e6bc 100644 --- a/gsw/__init__.py +++ b/gsw/__init__.py @@ -36,7 +36,7 @@ from .stability import * from .geostrophy import * from .utility import * -from .gibbs_function import gibbs, gibbs_ufunc +from .gibbs_function import gibbs from . import geostrophy from . import utility from . import stability diff --git a/gsw/gibbs_function.py b/gsw/gibbs_function.py index ca35a74..6457925 100644 --- a/gsw/gibbs_function.py +++ b/gsw/gibbs_function.py @@ -19,55 +19,46 @@ def _get_compiled_gibbs(): _gibbs_compiled = _get_compiled_gibbs() -def _set_gibbs(): - gibbs_compiled = _get_compiled_gibbs() - - def gibbs(ns, nt, np, SA, CT, p): - """ - Calculates the specific Gibbs free energy and derivatives up to order 2 - - Parameters - ---------- - ns : int - order of SA derivative - nt : int - order of CT derivative - np : int - order of p derivative - SA : array-like - Absolute Salinity, g/kg - CT : array-like - Conservative Temperature (ITS-90), degrees C - p : array-like - Sea pressure (absolute pressure minus 10.1325 dbar), dbar - - Returns - ------- - gibbs : specific Gibbs energy [J/kg] or its derivative - - - """ - if isinstance(SA, numpy.ndarray): - shape = SA.shape - assert CT.shape == SA.shape - assert p.shape == SA.shape - value = numpy.zeros(shape) - for k in numpy.ndindex(shape): - value[k] = gibbs_compiled(ns, nt, np, - c_double(SA[k]), c_double(CT[k]), c_double(p[k])) - else: - value = gibbs_compiled(ns, nt, np, - c_double(SA), c_double(CT), c_double(p)) - - return value - return gibbs - - def _basic_gibbs(ns, nt, np, SA, CT, p): """ this doc is lost during the ufunc-ification""" return _gibbs_compiled(ns, nt, np, c_double(SA), c_double(CT), c_double(p)) -gibbs = _set_gibbs() +_gibbs_ufunc = numpy.frompyfunc(_basic_gibbs, 6, 1) + + +def gibbs(ns, nt, np, SA, CT, p, **kwargs): + """Calculates the specific Gibbs free energy and derivatives up to order 2 + + Parameters + ---------- + ns : int + order of SA derivative + nt : int + order of CT derivative + np : int + order of p derivative + SA : array-like + Absolute Salinity, g/kg + CT : array-like + Conservative Temperature (ITS-90), degrees C + p : array-like + Sea pressure (absolute pressure minus 10.1325 dbar), dbar + + **kwargs: If either SA, CT or p is an nd.array then `gibbs` behaves + like a `ufunc` and **kwargs can be any of `ufunc` + keywords. See the :ref:`ufunc docs `. + + + Returns + ------- + gibbs : specific Gibbs energy [J/kg] or its derivative + + """ + if (isinstance(SA, numpy.ndarray) + or isinstance(CT, numpy.ndarray) + or isinstance(p, numpy.ndarray)): + return _gibbs_ufunc(ns, nt, np, SA, CT, p, **kwargs) -gibbs_ufunc = numpy.frompyfunc(_basic_gibbs, 6, 1) + return _gibbs_compiled(ns, nt, np, + c_double(SA), c_double(CT), c_double(p))