Validate the sto limits subbed into the OCPs give the correct voltage limits

CHANGELOG.md

@@ -2,6 +2,9 @@
 - Allow user-defined parameters to be added using the field ["Parameterisation"]["User-defined"] ([#44](https://github.com/pybamm-team/BPX/pull/44))
 - Added validation based on models: SPM, SPMe, DFN ([#34](https://github.com/pybamm-team/BPX/pull/34)). A warning will be produced if the user-defined model type does not match the parameter set (e.g., if the model is `SPM`, but the full DFN model parameters are provided).
 - Added support for well-mixed, blended electrodes that contain more than one active material ([#33](https://github.com/pybamm-team/BPX/pull/33))
+- Added validation of the STO limits subbed into the OCPs vs the upper/lower cut-off voltage limits for non-blended electrodes with the OCPs defined as functions ([#32](https://github.com/FaradayInstitution/BPX/pull/32)). The user can provide a tolerance by updating the settings variable `BPX.settings.tolerances["Voltage [V]"]` or by passing extra option `v_tol` to `parse_bpx_file()`, `parse_bpx_obj()` or `parse_bpx_str()` functions. Default value of the tolerance is 1 mV. The tolerance cannot be negative.
+- Added the target SOC check in `get_electrode_concentrations()` function. Raise a warning if the SOC is outside of [0,1] interval.
+- In `get_electrode_stoichiometries()` function, raise a warning instead of an error if the SOC is outside of [0,1] interval.
 - Added five parametrisation examples (two DFN parametrisation examples from About:Energy open-source release, blended electrode definition, user-defined 0th-order hysteresis, and SPM parametrisation).
 
 # [v0.3.1](https://github.com/pybamm-team/BPX/releases/tag/v0.3.1)

bpx/__init__.py

@@ -6,6 +6,7 @@
 from .interpolated_table import InterpolatedTable
 from .expression_parser import ExpressionParser
 from .function import Function
+from .validators import check_sto_limits
 from .schema import BPX
 from .parsers import parse_bpx_str, parse_bpx_obj, parse_bpx_file
 from .utilities import get_electrode_stoichiometries, get_electrode_concentrations

bpx/parsers.py

@@ -1,49 +1,70 @@
 from bpx import BPX
 
 
-def parse_bpx_file(filename: str) -> BPX:
+def parse_bpx_file(filename: str, v_tol: float = 0.001) -> BPX:
     """
     Parameters
     ----------
 
     filename: str
         a filepath to a bpx file
+    v_tol: float
+        absolute tolerance in [V] to validate the voltage limits, 1 mV by default
 
     Returns
     -------
     BPX:
         a parsed BPX model
     """
+    if v_tol < 0:
+        raise ValueError("v_tol should not be negative")
+
+    BPX.settings.tolerances["Voltage [V]"] = v_tol
+
     return BPX.parse_file(filename)
 
 
-def parse_bpx_obj(bpx: dict) -> BPX:
+def parse_bpx_obj(bpx: dict, v_tol: float = 0.001) -> BPX:
     """
     Parameters
     ----------
 
     bpx: dict
         a dict object in bpx format
+    v_tol: float
+        absolute tolerance in [V] to validate the voltage limits, 1 mV by default
 
     Returns
     -------
     BPX:
         a parsed BPX model
     """
+    if v_tol < 0:
+        raise ValueError("v_tol should not be negative")
+
+    BPX.settings.tolerances["Voltage [V]"] = v_tol
+
     return BPX.parse_obj(bpx)
 
 
-def parse_bpx_str(bpx: str) -> BPX:
+def parse_bpx_str(bpx: str, v_tol: float = 0.001) -> BPX:
     """
     Parameters
     ----------
 
     bpx: str
         a json formatted string in bpx format
+    v_tol: float
+        absolute tolerance in [V] to validate the voltage limits, 1 mV by default
 
     Returns
     -------
     BPX:
         a parsed BPX model
     """
+    if v_tol < 0:
+        raise ValueError("v_tol should not be negative")
+
+    BPX.settings.tolerances["Voltage [V]"] = v_tol
+
     return BPX.parse_raw(bpx)

bpx/schema.py

@@ -1,6 +1,6 @@
 from typing import List, Literal, Union, Dict, get_args
 from pydantic import BaseModel, Field, Extra, root_validator
-from bpx import Function, InterpolatedTable
+from bpx import Function, InterpolatedTable, check_sto_limits
 from warnings import warn
 
 FloatFunctionTable = Union[float, Function, InterpolatedTable]
@@ -10,6 +10,16 @@ class ExtraBaseModel(BaseModel):
     class Config:
         extra = Extra.forbid
 
+    class settings:
+        """
+        Class with BPX-related settings.
+        It might be worth moving it to a separate file if it grows bigger.
+        """
+
+        tolerances = {
+            "Voltage [V]": 1e-3,  # Absolute tolerance in [V] to validate the voltage limits
+        }
+
 
 class Header(ExtraBaseModel):
     bpx: float = Field(
@@ -191,7 +201,7 @@ class Particle(ExtraBaseModel):
     )
     maximum_concentration: float = Field(
         alias="Maximum concentration [mol.m-3]",
-        example=631040,
+        example=63104.0,
         description="Maximum concentration of lithium ions in particles",
     )
     particle_radius: float = Field(
@@ -338,6 +348,11 @@ class Parameterisation(ExtraBaseModel):
         alias="User-defined",
     )
 
+    # Reusable validators
+    _sto_limit_validation = root_validator(skip_on_failure=True, allow_reuse=True)(
+        check_sto_limits
+    )
+
 
 class ParameterisationSPM(ExtraBaseModel):
     cell: Cell = Field(
@@ -354,6 +369,11 @@ class ParameterisationSPM(ExtraBaseModel):
         alias="User-defined",
     )
 
+    # Reusable validators
+    _sto_limit_validation = root_validator(skip_on_failure=True, allow_reuse=True)(
+        check_sto_limits
+    )
+
 
 class BPX(ExtraBaseModel):
     header: Header = Field(

bpx/utilities.py

@@ -1,3 +1,6 @@
+from warnings import warn
+
+
 def get_electrode_stoichiometries(target_soc, bpx):
     """
     Calculate individual electrode stoichiometries at a particular target
@@ -16,7 +19,7 @@ def get_electrode_stoichiometries(target_soc, bpx):
         The electrode stoichiometries that give the target state of charge
     """
     if target_soc < 0 or target_soc > 1:
-        raise ValueError("Target SOC should be between 0 and 1")
+        warn("Target SOC should be between 0 and 1")
 
     sto_n_min = bpx.parameterisation.negative_electrode.minimum_stoichiometry
     sto_n_max = bpx.parameterisation.negative_electrode.maximum_stoichiometry
@@ -47,6 +50,9 @@ def get_electrode_concentrations(target_soc, bpx):
     c_n, c_p
         The electrode concentrations that give the target state of charge
     """
+    if target_soc < 0 or target_soc > 1:
+        warn("Target SOC should be between 0 and 1")
+
     c_n_max = bpx.parameterisation.negative_electrode.maximum_concentration
     c_p_max = bpx.parameterisation.positive_electrode.maximum_concentration
 

bpx/validators.py

@@ -0,0 +1,47 @@
+from warnings import warn
+
+
+def check_sto_limits(cls, values):
+    """
+    Validates that the STO limits subbed into the OCPs give the correct voltage limits.
+    Works if both OCPs are defined as functions.
+    Blended electrodes are not supported.
+    This is a reusable validator to be used for both DFN/SPMe and SPM parameter sets.
+    """
+
+    try:
+        ocp_n = values.get("negative_electrode").ocp.to_python_function()
+        ocp_p = values.get("positive_electrode").ocp.to_python_function()
+    except AttributeError:
+        # OCPs defined as interpolated tables or one of the electrodes is blended; do nothing
+        return values
+
+    sto_n_min = values.get("negative_electrode").minimum_stoichiometry
+    sto_n_max = values.get("negative_electrode").maximum_stoichiometry
+    sto_p_min = values.get("positive_electrode").minimum_stoichiometry
+    sto_p_max = values.get("positive_electrode").maximum_stoichiometry
+    V_min = values.get("cell").lower_voltage_cutoff
+    V_max = values.get("cell").upper_voltage_cutoff
+
+    # Voltage tolerance from `settings` data class
+    tol = cls.settings.tolerances["Voltage [V]"]
+
+    # Checks the maximum voltage estimated from STO
+    V_max_sto = ocp_p(sto_p_min) - ocp_n(sto_n_max)
+    if V_max_sto - V_max > tol:
+        warn(
+            f"The maximum voltage computed from the STO limits ({V_max_sto} V) "
+            f"is higher than the upper voltage cut-off ({V_max} V) "
+            f"with the absolute tolerance v_tol = {tol} V"
+        )
+
+    # Checks the minimum voltage estimated from STO
+    V_min_sto = ocp_p(sto_p_max) - ocp_n(sto_n_min)
+    if V_min_sto - V_min < -tol:
+        warn(
+            f"The minimum voltage computed from the STO limits ({V_min_sto} V) "
+            f"is less than the lower voltage cut-off ({V_min} V) "
+            f"with the absolute tolerance v_tol = {tol} V"
+        )
+
+    return values

tests/test_parsers.py

@@ -0,0 +1,123 @@
+import unittest
+import warnings
+import copy
+
+from bpx import BPX, parse_bpx_file, parse_bpx_obj, parse_bpx_str
+
+
+class TestParsers(unittest.TestCase):
+    def setUp(self):
+        base = """
+            {
+                "Header": {
+                        "BPX": 0.1,
+                        "Title": "Parameterisation example of an NMC111|graphite 12.5 Ah pouch cell",
+                        "Model": "DFN"
+                },
+                "Parameterisation": {
+                        "Cell": {
+                            "Ambient temperature [K]": 298.15,
+                            "Initial temperature [K]": 298.15,
+                            "Reference temperature [K]": 298.15,
+                            "Lower voltage cut-off [V]": 2.7,
+                            "Upper voltage cut-off [V]": 4.2,
+                            "Nominal cell capacity [A.h]": 12.5,
+                            "Specific heat capacity [J.K-1.kg-1]": 913,
+                            "Thermal conductivity [W.m-1.K-1]": 2.04,
+                            "Density [kg.m-3]": 1847,
+                            "Electrode area [m2]": 0.016808,
+                            "Number of electrode pairs connected in parallel to make a cell": 34,
+                            "External surface area [m2]": 0.0379,
+                            "Volume [m3]": 0.000128
+                        },
+                        "Electrolyte": {
+                            "Initial concentration [mol.m-3]": 1000,
+                            "Cation transference number": 0.2594,
+                            "Conductivity [S.m-1]": "0.1297 * (x / 1000) ** 3 - 2.51 * (x / 1000) ** 1.5 + 3.329 * (x / 1000)",
+                            "Diffusivity [m2.s-1]": "8.794e-11 * (x / 1000) ** 2 - 3.972e-10 * (x / 1000) + 4.862e-10",
+                            "Conductivity activation energy [J.mol-1]": 17100,
+                            "Diffusivity activation energy [J.mol-1]": 17100
+                        },
+                        "Negative electrode": {
+                            "Particle radius [m]": 4.12e-06,
+                            "Thickness [m]": 5.62e-05,
+                            "Diffusivity [m2.s-1]": 2.728e-14,
+                            "OCP [V]":
+                                "9.47057878e-01 * exp(-1.59418743e+02  * x) - 3.50928033e+04 +
+                                1.64230269e-01 * tanh(-4.55509094e+01 * (x - 3.24116012e-02 )) +
+                                3.69968491e-02 * tanh(-1.96718868e+01 * (x - 1.68334476e-01)) +
+                                1.91517003e+04 * tanh(3.19648312e+00 * (x - 1.85139824e+00)) +
+                                5.42448511e+04 * tanh(-3.19009848e+00 * (x - 2.01660395e+00))",
+                            "Entropic change coefficient [V.K-1]":
+                                "(-0.1112 * x + 0.02914 + 0.3561 * exp(-((x - 0.08309) ** 2) / 0.004616)) / 1000",
+                            "Conductivity [S.m-1]": 0.222,
+                            "Surface area per unit volume [m-1]": 499522,
+                            "Porosity": 0.253991,
+                            "Transport efficiency": 0.128,
+                            "Reaction rate constant [mol.m-2.s-1]": 5.199e-06,
+                            "Minimum stoichiometry": 0.005504,
+                            "Maximum stoichiometry": 0.75668,
+                            "Maximum concentration [mol.m-3]": 29730,
+                            "Diffusivity activation energy [J.mol-1]": 30000,
+                            "Reaction rate constant activation energy [J.mol-1]": 55000
+                        },
+                        "Positive electrode": {
+                            "Particle radius [m]": 4.6e-06,
+                            "Thickness [m]": 5.23e-05,
+                            "Diffusivity [m2.s-1]": 3.2e-14,
+                            "OCP [V]":
+                                "-3.04420906 * x + 10.04892207 - 0.65637536 * tanh(-4.02134095 * (x - 0.80063948)) +
+                                4.24678547 * tanh(12.17805062 * (x - 7.57659337)) - 0.3757068 * tanh(59.33067782 * (x - 0.99784492))",
+                            "Entropic change coefficient [V.K-1]": -1e-4,
+                            "Conductivity [S.m-1]": 0.789,
+                            "Surface area per unit volume [m-1]": 432072,
+                            "Porosity": 0.277493,
+                            "Transport efficiency": 0.1462,
+                            "Reaction rate constant [mol.m-2.s-1]": 2.305e-05,
+                            "Minimum stoichiometry": 0.42424,
+                            "Maximum stoichiometry": 0.96210,
+                            "Maximum concentration [mol.m-3]": 46200,
+                            "Diffusivity activation energy [J.mol-1]": 15000,
+                            "Reaction rate constant activation energy [J.mol-1]": 35000
+                        },
+                        "Separator": {
+                            "Thickness [m]": 2e-05,
+                            "Porosity": 0.47,
+                            "Transport efficiency": 0.3222
+                        }
+                }
+            }
+            """
+        self.base = base.replace("\n", "")
+
+    def test_negative_v_tol_file(self):
+        with self.assertRaisesRegex(
+            ValueError,
+            "v_tol should not be negative",
+        ):
+            parse_bpx_file("filename", v_tol=-0.001)
+
+    def test_negative_v_tol_object(self):
+        bpx_obj = {"BPX": 1.0}
+        with self.assertRaisesRegex(
+            ValueError,
+            "v_tol should not be negative",
+        ):
+            parse_bpx_obj(bpx_obj, v_tol=-0.001)
+
+    def test_negative_v_tol_string(self):
+        with self.assertRaisesRegex(
+            ValueError,
+            "v_tol should not be negative",
+        ):
+            parse_bpx_str("String", v_tol=-0.001)
+
+    def test_parse_string(self):
+        test = copy.copy(self.base)
+        with self.assertWarns(UserWarning):
+            parse_bpx_str(test)
+
+    def test_parse_string_tolerance(self):
+        warnings.filterwarnings("error")  # Treat warnings as errors
+        test = copy.copy(self.base)
+        parse_bpx_str(test, v_tol=0.002)