Skip to content

Commit

Permalink
Fix issues with opt_k_carbonic 7-8
Browse files Browse the repository at this point in the history
mvdh7 committed Nov 17, 2023

Verified

This commit was created on GitHub.com and signed with GitHub’s verified signature.
1 parent bb0d553 commit 2c27953
Showing 4 changed files with 68 additions and 3 deletions.
7 changes: 5 additions & 2 deletions PyCO2SYS/equilibria/p1atm.py
Original file line number Diff line number Diff line change
@@ -400,7 +400,6 @@ def kH2CO3_SWS_HM_DM87(TempK, Sal):
return K1, K2


@np.errstate(divide="ignore", invalid="ignore") # because Sal=0 gives log10(Sal)=-inf
def kH2CO3_NBS_MCHP73(TempK, Sal):
"""Carbonic acid dissociation constants following MCHP73."""
# === CO2SYS.m comments: =======
@@ -409,6 +408,8 @@ def kH2CO3_NBS_MCHP73(TempK, Sal):
# I.e., these are the original Mehrbach dissociation constants.
# The 2s precision in pK1 is .005, or 1.2% in K1.
# The 2s precision in pK2 is .008, or 2% in K2.
Sal = np.where(Sal < 1e-16, 1e-16, Sal)
# ^ added in v1.8.3, because Sal=0 gives log10(Sal)=-inf
pK1 = (
-13.7201
+ 0.031334 * TempK
@@ -427,7 +428,9 @@ def kH2CO3_NBS_MCHP73(TempK, Sal):
- 8.0944e-4 * Sal * TempK
- 5617.11 * np.log10(Sal) / TempK
+ 2.136 * Sal / TempK
) # pK2 is not defined for Sal=0, since log10(0)=-inf
) # pK2 is not defined for Sal=0, since log10(0)=-inf, but since v1.8.3 we
# return the value for Sal=1e-16 instead (this option shouldn't be used in
# such low salinities anyway, it's only valid above 19!)
K2 = 10.0**-pK2 # this is on the NBS scale
return K1, K2

4 changes: 3 additions & 1 deletion PyCO2SYS/solve/__init__.py
Original file line number Diff line number Diff line change
@@ -560,7 +560,9 @@ def get_lnpCO2(
pressure_atmosphere=pressure_atmosphere,
opt_pressured_kCO2=opt_pressured_kCO2,
)
fCO2 = get.fCO2fromTATC(alkalinity, dic, totals, k_constants)
fCO2 = get.fCO2fromTATC(
alkalinity - totals["PengCorrection"], dic, totals, k_constants
)
pCO2 = convert.fCO2_to_pCO2(fCO2, k_constants)
return np.log(pCO2)

1 change: 1 addition & 0 deletions docs/versions.md
Original file line number Diff line number Diff line change
@@ -46,6 +46,7 @@ Adds atmospheric pressure input for *p*CO<sub>2</sub>-*f*CO<sub>2</sub>-*x*CO<su
***Technical***

* Updated from building with setup.py to pyproject.toml.
* Updated `pyco2.equilibria.p1atm.kH2CO3_NBS_MCHP73` (used for `opt_k_carbonic` options `6` and `7`) to update any salinity values less than 10<sup>–16</sup> to be 10<sup>–16</sup>, because zero salinities give a NaN for <i>K</i><sub>2</sub>, which causes autograd problems. This should not make any practical difference, because the parameterisation is only valid for salinities above 19.

### 1.8.2 (19 January 2023)

59 changes: 59 additions & 0 deletions tests/test_dlnpCO2_dT.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,59 @@
import PyCO2SYS as pyco2
import numpy as np

# Simulate arguments
rng = np.random.default_rng(1)
npts = 1000
alkalinity = rng.normal(loc=2300, scale=50, size=npts)
dic = rng.normal(loc=2150, scale=50, size=npts)
salinity = rng.uniform(low=0, high=40, size=npts)
temperature = rng.uniform(low=-2, high=35, size=npts)
pressure = rng.uniform(low=0, high=10000, size=npts)
temperature_out = rng.uniform(low=-2, high=35, size=npts)
pressure_out = rng.uniform(low=0, high=10000, size=npts)
total_silicate = rng.uniform(low=0, high=100, size=npts)
total_phosphate = rng.uniform(low=0, high=100, size=npts)
total_ammonia = rng.uniform(low=0, high=100, size=npts)
total_sulfide = rng.uniform(low=0, high=100, size=npts)
opt_k_carbonic = rng.integers(low=1, high=18, size=npts)

# Solve
results = pyco2.sys(
par1=alkalinity,
par2=dic,
par1_type=1,
par2_type=2,
salinity=salinity,
temperature=temperature,
pressure=pressure,
temperature_out=temperature_out,
pressure_out=pressure_out,
total_silicate=total_silicate,
total_phosphate=total_phosphate,
total_ammonia=total_ammonia,
total_sulfide=total_sulfide,
opt_k_carbonic=opt_k_carbonic,
opt_pressured_kCO2=0,
grads_of=["pCO2"],
grads_wrt=["temperature"],
)

# Extract values
dlnpCO2_dT_autograd = results["dlnpCO2_dT"]
dlnpCO2_dT_ffd = results["d_pCO2__d_temperature"] / results["pCO2"]
diff = dlnpCO2_dT_autograd - dlnpCO2_dT_ffd
abs_diff = np.abs(diff)
max_abs_diff = np.max(abs_diff)
diff_pct = 200 * diff / (dlnpCO2_dT_autograd + dlnpCO2_dT_ffd)
abs_diff_pct = np.abs(diff_pct)
max_abs_diff_pct = np.max(abs_diff_pct)


def test_dlnpCO2_dT():
"""Are all of the autograd values the same as the forward finite difference values
to within 0.02%?
"""
assert np.allclose(dlnpCO2_dT_autograd, dlnpCO2_dT_ffd, rtol=0.02 / 100, atol=0)


test_dlnpCO2_dT()

0 comments on commit 2c27953

Please sign in to comment.