From 51c106611e7c82e542bace773f1936578865ec3c Mon Sep 17 00:00:00 2001 From: Matthew Humphreys Date: Fri, 13 Dec 2024 22:55:25 +0100 Subject: [PATCH] v3.1.1 tests working --- PyCO2SYS/engine/system.py | 12 +- PyCO2SYS/equilibria/p1atm.py | 44 +++- .../compare_equilibrium_constants_v3_1_1.py | 235 +++++++++++++----- 3 files changed, 220 insertions(+), 71 deletions(-) diff --git a/PyCO2SYS/engine/system.py b/PyCO2SYS/engine/system.py index 81fb12b1..fe5c2be2 100644 --- a/PyCO2SYS/engine/system.py +++ b/PyCO2SYS/engine/system.py @@ -453,7 +453,13 @@ ), ), 17: dict( - k_H2CO3_sws_1atm=equilibria.p1atm.k_H2CO3_sws_WMW14, + # k_H2CO3_sws_1atm=equilibria.p1atm.k_H2CO3_sws_WMW14, + # ^ although the above should work, it gives slightly different answers than + # the conversion below, and below is consistent with the MATLAB implementation + k_H2CO3_total_1atm=equilibria.p1atm.k_H2CO3_total_WMW14, + k_H2CO3_sws_1atm=lambda k_H2CO3_total_1atm, tot_to_sws_1atm: ( + k_H2CO3_total_1atm * tot_to_sws_1atm + ), k_HCO3_total_1atm=equilibria.p1atm.k_HCO3_total_SB21, k_HCO3_sws_1atm=lambda k_HCO3_total_1atm, tot_to_sws_1atm: ( k_HCO3_total_1atm * tot_to_sws_1atm @@ -520,13 +526,13 @@ 3: dict(k_HSO4_free_1atm=equilibria.p1atm.k_HSO4_free_WM13), } get_funcs_opts["opt_k_NH3"] = { - 1: dict(k_NH3_sws_1atm=equilibria.p1atm.k_NH3_sws_YM95), - 2: dict( + 1: dict( k_NH3_total_1atm=equilibria.p1atm.k_NH3_total_CW95, k_NH3_sws_1atm=lambda k_NH3_total_1atm, tot_to_sws_1atm: ( k_NH3_total_1atm * tot_to_sws_1atm ), ), + 2: dict(k_NH3_sws_1atm=equilibria.p1atm.k_NH3_sws_YM95), } get_funcs_opts["opt_k_Si"] = { 1: dict(k_Si_sws_1atm=equilibria.p1atm.k_Si_sws_YM95), diff --git a/PyCO2SYS/equilibria/p1atm.py b/PyCO2SYS/equilibria/p1atm.py index b7a7ff0e..da624f74 100644 --- a/PyCO2SYS/equilibria/p1atm.py +++ b/PyCO2SYS/equilibria/p1atm.py @@ -1543,21 +1543,57 @@ def k_HCO3_sws_WMW14(temperature, salinity): return 10.0**-pK2 -def k_H2CO3_TOT_WMW14(TempK, salinity): - """Carbonic acid dissociation constants, Total scale, following WM13/WMW14.""" +def k_H2CO3_total_WMW14(temperature, salinity): + """First carbonic acid dissociation constant following WM13/WMW14. + Used when opt_k_carbonic = 17. + + Parameters + ---------- + temperature : float + Temperature in °C. + salinity : float + Practical salinity. + + Returns + ------- + float + HCO3 dissociation constant. + """ + TempK = convert.celsius_to_kelvin(temperature) # Coefficients from the corrigendum document [WMW14] - pK10, pK20 = _kH2CO3_WMW14(TempK) + pK10 = _kH2CO3_WMW14(TempK)[0] A1 = 13.568513 * salinity**0.5 + 0.031645 * salinity - 5.3834e-5 * salinity**2 B1 = -539.2304 * salinity**0.5 - 5.635 * salinity C1 = -2.0901396 * salinity**0.5 pK1 = pK10 + A1 + B1 / TempK + C1 * np.log(TempK) K1 = 10.0**-pK1 + return K1 + + +def k_HCO3_total_WMW14(temperature, salinity): + """Second carbonic acid dissociation constant following WM13/WMW14. + + Parameters + ---------- + temperature : float + Temperature in °C. + salinity : float + Practical salinity. + + Returns + ------- + float + HCO3 dissociation constant. + """ + TempK = convert.celsius_to_kelvin(temperature) + # Coefficients from the corrigendum document [WMW14] + pK20 = _kH2CO3_WMW14(TempK)[1] A2 = 21.389248 * salinity**0.5 + 0.12452358 * salinity - 3.7447e-4 * salinity**2 B2 = -787.3736 * salinity**0.5 - 19.84233 * salinity C2 = -3.3773006 * salinity**0.5 pK2 = pK20 + A2 + B2 / TempK + C2 * np.log(TempK) K2 = 10.0**-pK2 - return K1, K2 + return K2 def k_H2CO3_FREE_WMW14(TempK, salinity): diff --git a/manuscript/compare_equilibrium_constants_v3_1_1.py b/manuscript/compare_equilibrium_constants_v3_1_1.py index 2f83ad24..6e7a1427 100644 --- a/manuscript/compare_equilibrium_constants_v3_1_1.py +++ b/manuscript/compare_equilibrium_constants_v3_1_1.py @@ -1,84 +1,191 @@ import warnings -import pandas as pd, numpy as np -import PyCO2SYS as pyco2 + +import numpy as np +import pandas as pd + +from PyCO2SYS import CO2System # Import MATLAB results and recalculate with PyCO2SYS matlab = pd.read_csv("manuscript/results/compare_equilibrium_constants_v3_1_1.csv") -python = pd.DataFrame( - pyco2.sys( - matlab.PAR1.values, - matlab.PAR2.values, - matlab.PAR1TYPE.values, - matlab.PAR2TYPE.values, - salinity=matlab.SAL.values, - temperature=matlab.TEMPIN.values, - temperature_out=matlab.TEMPOUT.values, - pressure=matlab.PRESIN.values, - pressure_out=matlab.PRESOUT.values, - opt_pH_scale=matlab.pHSCALEIN.values, - opt_gas_constant=3, - opt_k_carbonic=matlab.K1K2CONSTANTS.values, - opt_total_borate=matlab.BORON.values, - opt_k_bisulfate=matlab.KSO4CONSTANT.values, - opt_k_fluoride=matlab.KFCONSTANT.values, - total_phosphate=matlab.PO4.values, - total_silicate=matlab.SI.values, - total_ammonia=matlab.TNH4.values, - total_sulfide=matlab.TH2S.values, - opt_buffers_mode=0, - ) -) def test_equilibrium_constants(): - for m, p in ( - ("K0input", "k_CO2"), - ("K0output", "k_CO2_out"), - ("K1input", "k_carbonic_1"), - ("K1output", "k_carbonic_1_out"), - ("K2input", "k_carbonic_2"), - ("K2output", "k_carbonic_2_out"), - ("KWinput", "k_water"), - ("KWoutput", "k_water_out"), - ("KP1input", "k_phosphoric_1"), - ("KP1output", "k_phosphoric_1_out"), - ("KP2input", "k_phosphoric_2"), - ("KP2output", "k_phosphoric_2_out"), - ("KP3input", "k_phosphoric_3"), - ("KP3output", "k_phosphoric_3_out"), - ("KFinput", "k_fluoride"), - ("KFoutput", "k_fluoride_out"), - ("KSinput", "k_bisulfate"), - ("KSoutput", "k_bisulfate_out"), - ("KSiinput", "k_silicate"), - ("KSioutput", "k_silicate_out"), - ("KBinput", "k_borate"), - ("KBoutput", "k_borate_out"), - ("KNH4input", "k_ammonia"), - ("KNH4output", "k_ammonia_out"), - ("KH2Sinput", "k_sulfide"), - ("KH2Soutput", "k_sulfide_out"), + m_to_p = ( + ("K0", "k_CO2_1atm"), + ("K1", "k_H2CO3"), + ("K2", "k_HCO3"), + ("KW", "k_H2O"), + ("KP1", "k_H3PO4"), + ("KP2", "k_H2PO4"), + ("KP3", "k_HPO4"), + ("KF", "k_HF_free"), + ("KS", "k_HSO4_free"), + ("KSi", "k_Si"), + ("KB", "k_BOH3"), + ("KNH4", "k_NH3"), + ("KH2S", "k_H2S"), + ) + svars = [p for m, p in m_to_p] + for g, group in matlab.groupby( + ["K1K2CONSTANTS", "pHSCALEIN", "KSO4CONSTANT", "BORON", "KFCONSTANT"] ): - with warnings.catch_warnings(): - warnings.filterwarnings("ignore", category=RuntimeWarning) - pk_matlab = np.where( - matlab[m].values == 0, -999.9, -np.log10(matlab[m].values) + # Set up arguments dicts + values_in = dict( + temperature=group.TEMPIN.values, + pressure=group.PRESIN.values, + salinity=group.SAL.values, + ) + values_out = dict( + temperature=group.TEMPOUT.values, + pressure=group.PRESOUT.values, + salinity=group.SAL.values, + ) + opts = dict( + opt_k_carbonic=g[0], + opt_pH_scale=g[1], + opt_k_HSO4=g[2], + opt_total_borate=g[3], + opt_k_HF=g[4], + opt_gas_constant=3, + ) + # Deal with GEOSECS and freshwater weirdness + if g[0] == 6: + opts.update( + dict( + opt_k_BOH3=2, + opt_factor_k_BOH3=2, + opt_factor_k_H2CO3=3, + opt_factor_k_HCO3=3, + ) + ) + elif g[0] == 7: + opts.update( + dict( + opt_fH=2, + opt_k_BOH3=2, + opt_k_H2O=2, + opt_k_phosphate=2, + opt_k_Si=2, + opt_factor_k_BOH3=2, + opt_factor_k_H2CO3=3, + opt_factor_k_HCO3=3, + ) + ) + elif g[0] == 8: + values_in.update(dict(salinity=0.0)) + values_out.update(dict(salinity=0.0)) + opts.update( + dict( + opt_fH=3, + opt_k_H2O=3, + opt_factor_k_H2O=2, + opt_factor_k_H2CO3=2, + opt_factor_k_HCO3=2, + ) ) - pk_python = np.where( - python[p].values == 0, -999.9, -np.log10(python[p].values) + # Solve under input and output conditions + sys_in = CO2System(values=values_in, opts=opts) + sys_in.solve(svars) + sys_out = CO2System(values=values_out, opts=opts) + sys_out.solve(svars) + # Compare MATLAB with Python + for m, p in m_to_p: + with warnings.catch_warnings(): + warnings.filterwarnings("ignore", category=RuntimeWarning) + pk_matlab_in = np.where( + group[m + "input"].values == 0, + -999.9, + -np.log10(group[m + "input"].values), + ) + pk_python_in = np.where( + sys_in.values[p] == 0, -999.9, -np.log10(sys_in.values[p]) + ) + pk_matlab_out = np.where( + group[m + "output"].values == 0, + -999.9, + -np.log10(group[m + "output"].values), + ) + pk_python_out = np.where( + sys_out.values[p] == 0, -999.9, -np.log10(sys_out.values[p]) + ) + # These terms are not included when opt_k_carbonic == 6 + if g[0] == 6 and p in [ + "k_H2O", + "k_H3PO4", + "k_H2PO4", + "k_HPO4", + "k_Si", + "k_NH3", + "k_H2S", + ]: + pk_python_in[:] = -999.9 + pk_python_out[:] = -999.9 + # These terms are not included when opt_k_carbonic == 7 + if g[0] == 7 and p in [ + "k_NH3", + "k_H2S", + ]: + pk_python_in[:] = -999.9 + pk_python_out[:] = -999.9 + # These terms are not included when opt_k_carbonic == 8 + if g[0] == 8 and p in [ + "k_H3PO4", + "k_H2PO4", + "k_HPO4", + "k_Si", + "k_BOH3", + "k_NH3", + "k_H2S", + ]: + pk_python_in[:] = -999.9 + pk_python_out[:] = -999.9 + assert np.all( + np.isclose( + pk_matlab_in, + pk_python_in, + rtol=1e-12, + atol=1e-16, + ) + ) + assert np.all( + np.isclose( + pk_matlab_out, + pk_python_out, + rtol=1e-12, + atol=1e-16, + ) ) - assert np.all(np.isclose(pk_matlab, pk_python, rtol=1e-12, atol=1e-16)) def test_total_salts(): - for m, p in ( + m_to_p = ( ("TS", "total_sulfate"), ("TF", "total_fluoride"), ("TB", "total_borate"), - ): - assert np.all( - np.isclose(matlab[m].values, python[p].values, rtol=1e-12, atol=1e-16) + ) + svars = [p for m, p in m_to_p] + for g, group in matlab.groupby(["K1K2CONSTANTS", "BORON"]): + # Set up arguments dicts + values = dict(salinity=group.SAL.values) + opts = dict( + opt_k_carbonic=g[0], + opt_total_borate=g[1], ) + # Deal with GEOSECS and freshwater weirdness + if g[0] == 6: + opts.update(dict(opt_total_borate=4)) + elif g[0] == 7: + opts.update(dict(opt_total_borate=4)) + # Solve + sys = CO2System(values=values, opts=opts) + sys.solve(svars) + # Compare MATLAB with Python + for m, p in m_to_p: + python = sys.values[p] + # These terms are not included when opt_k_carbonic == 8 + if g[0] == 8 and p in ["total_sulfate", "total_fluoride", "total_borate"]: + python[:] = 0.0 + assert np.all(np.isclose(group[m].values, python, rtol=1e-12, atol=1e-16)) # test_equilibrium_constants()