Skip to content

Commit

Permalink
v3.1.1 tests working
Browse files Browse the repository at this point in the history
  • Loading branch information
@mvdh7
mvdh7 committed Dec 13, 2024
1 parent 63d1536 commit 51c1066
Show file tree
Hide file tree
Showing 3 changed files with 220 additions and 71 deletions.
12 changes: 9 additions & 3 deletions PyCO2SYS/engine/system.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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),
Expand Down
44 changes: 40 additions & 4 deletions PyCO2SYS/equilibria/p1atm.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand Down
235 changes: 171 additions & 64 deletions manuscript/compare_equilibrium_constants_v3_1_1.py
View file
Original file line number Diff line number Diff line change
@@ -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()
Expand Down

0 comments on commit 51c1066

Please sign in to comment.