From f3faca473c419f41e6cd7b1fdec3a89a2f094136 Mon Sep 17 00:00:00 2001 From: Caitlin Parke Date: Fri, 9 Aug 2024 16:33:10 -0400 Subject: [PATCH 01/30] Added ECM + split OCV model --- .../lithium_ion/__init__.py | 3 +- .../lithium_ion/basic_ecm_split_OCV.py | 87 +++++++++++++++++++ 2 files changed, 89 insertions(+), 1 deletion(-) create mode 100644 src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py diff --git a/src/pybamm/models/full_battery_models/lithium_ion/__init__.py b/src/pybamm/models/full_battery_models/lithium_ion/__init__.py index b02868dbe9..f98da10b73 100644 --- a/src/pybamm/models/full_battery_models/lithium_ion/__init__.py +++ b/src/pybamm/models/full_battery_models/lithium_ion/__init__.py @@ -24,8 +24,9 @@ from .Yang2017 import Yang2017 from .mpm import MPM from .msmr import MSMR +from .basic_ecm_split_OCV import ECMsplitOCV __all__ = ['Yang2017', 'base_lithium_ion_model', 'basic_dfn', 'basic_dfn_composite', 'basic_dfn_half_cell', 'basic_spm', 'dfn', 'electrode_soh', 'electrode_soh_half_cell', 'mpm', 'msmr', - 'newman_tobias', 'spm', 'spme'] + 'newman_tobias', 'spm', 'spme', 'ecm_split_ocv'] diff --git a/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py b/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py new file mode 100644 index 0000000000..5ad9574ab2 --- /dev/null +++ b/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py @@ -0,0 +1,87 @@ +# +# Equivalent Circuit Model with split OCV +# +import pybamm +from .base_lithium_ion_model import BaseModel + + +class ECMsplitOCV(BaseModel): + """Basic Equivalent Circuit Model that uses two OCV functions + for each electrode from the OCV function from Lithium ion parameter sets. + This class differs from the :class: pybamm.equivalent_circuit.Thevenin() due + to dual OCV functions to make up the voltage from each electrode. + + Parameters + ---------- + name: str, optional + The name of the model. + """ + + def __init__(self, name="ECM with split OCV"): + super().__init__({}, name) + # TODO citations + param = self.param + + ###################### + # Variables + ###################### + # All variables are only time-dependent + # No domain definition needed + + c_n = pybamm.Variable("Negative particle SOC") + c_p = pybamm.Variable("Positive particle SOC") + Q = pybamm.Variable("Discharge capacity [A.h]") + V = pybamm.Variable("Voltage [V]") + + # model is isothermal + T = param.T_init + I = param.current_with_time + + # Capacity equation + self.rhs[Q] = I / 3600 + self.initial_conditions[Q] = pybamm.Scalar(0) + + # Capacity in each electrode + # TODO specify capcity for negative and positive electrodes + # may be user-defined + q_n = 1 # Ah + q_p = 1 # Ah + Qn = q_n + Qp = q_p + + # State of charge electrode equations + self.rhs[c_n] = - I / Qn / 3600 + self.rhs[c_p] = I / Qp / 3600 + self.initial_conditions[c_n] = param.n.prim.c_init_av / param.n.prim.c_max + self.initial_conditions[c_p] = param.p.prim.c_init_av / param.p.prim.c_max + + # OCV's for the electrodes + Un = param.n.prim.U(c_n, T) + Up = param.p.prim.U(c_p, T) + + # IR resistance, hard coded for now + IR = 0.1 + V = Up - Un - IR + + self.variables = { + "Negative particle SOC": c_n, + "Positive particle SOC": c_p, + "Current [A]": I, + "Discharge capacity [A.h]": Q, + "Voltage [V]": V, + "Times [s]": pybamm.t, + "Positive electrode potential [V]": Up, + "Negative electrode potential [V]": Un, + "Current variable [A]": I, + "Current function [A]": I, + } + + # Events specify points at which a solution should terminate + self.events += [ + pybamm.Event("Minimum voltage [V]", V - param.voltage_low_cut), + pybamm.Event("Maximum voltage [V]", param.voltage_high_cut - V), + pybamm.Event("Maximum Negative Electrode SOC", 0.999 - c_n), + pybamm.Event("Maximum Positive Electrode SOC", 0.999 - c_p), + pybamm.Event("Minimum Negative Electrode SOC", c_n - 0.0001), + pybamm.Event("Minimum Positive Electrode SOC", c_p - 0.0001), + ] From f93e67b2fcaf2a183cb7fd63774196b698807a0f Mon Sep 17 00:00:00 2001 From: Caitlin Parke Date: Fri, 9 Aug 2024 16:39:33 -0400 Subject: [PATCH 02/30] pre-commit changes --- .../full_battery_models/lithium_ion/basic_ecm_split_OCV.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py b/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py index 5ad9574ab2..7d8d933a20 100644 --- a/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py +++ b/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py @@ -44,13 +44,13 @@ def __init__(self, name="ECM with split OCV"): # Capacity in each electrode # TODO specify capcity for negative and positive electrodes # may be user-defined - q_n = 1 # Ah - q_p = 1 # Ah + q_n = 1 # Ah + q_p = 1 # Ah Qn = q_n Qp = q_p # State of charge electrode equations - self.rhs[c_n] = - I / Qn / 3600 + self.rhs[c_n] = -I / Qn / 3600 self.rhs[c_p] = I / Qp / 3600 self.initial_conditions[c_n] = param.n.prim.c_init_av / param.n.prim.c_max self.initial_conditions[c_p] = param.p.prim.c_init_av / param.p.prim.c_max From 83239a3d48b9ad5ce1564f16f299d73f590389e2 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 9 Aug 2024 20:43:34 +0000 Subject: [PATCH 03/30] style: pre-commit fixes --- .../full_battery_models/lithium_ion/basic_ecm_split_OCV.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py b/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py index 5ad9574ab2..7d8d933a20 100644 --- a/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py +++ b/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py @@ -44,13 +44,13 @@ def __init__(self, name="ECM with split OCV"): # Capacity in each electrode # TODO specify capcity for negative and positive electrodes # may be user-defined - q_n = 1 # Ah - q_p = 1 # Ah + q_n = 1 # Ah + q_p = 1 # Ah Qn = q_n Qp = q_p # State of charge electrode equations - self.rhs[c_n] = - I / Qn / 3600 + self.rhs[c_n] = -I / Qn / 3600 self.rhs[c_p] = I / Qp / 3600 self.initial_conditions[c_n] = param.n.prim.c_init_av / param.n.prim.c_max self.initial_conditions[c_p] = param.p.prim.c_init_av / param.p.prim.c_max From 2a945f280c428610fd74c008eb9bee97baf33f9f Mon Sep 17 00:00:00 2001 From: Caitlin Parke Date: Mon, 12 Aug 2024 11:17:11 -0400 Subject: [PATCH 04/30] Added parameter objects in model --- .../lithium_ion/basic_ecm_split_OCV.py | 20 +++++++++---------- 1 file changed, 9 insertions(+), 11 deletions(-) diff --git a/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py b/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py index 7d8d933a20..4cb0730958 100644 --- a/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py +++ b/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py @@ -42,16 +42,12 @@ def __init__(self, name="ECM with split OCV"): self.initial_conditions[Q] = pybamm.Scalar(0) # Capacity in each electrode - # TODO specify capcity for negative and positive electrodes - # may be user-defined - q_n = 1 # Ah - q_p = 1 # Ah - Qn = q_n - Qp = q_p + Q_n = pybamm.Parameter("Negative electrode capacity [A.h]") + Q_p = pybamm.Parameter("Positive electrode capacity [A.h]") # State of charge electrode equations - self.rhs[c_n] = -I / Qn / 3600 - self.rhs[c_p] = I / Qp / 3600 + self.rhs[c_n] = -I / Q_n / 3600 + self.rhs[c_p] = I / Q_p / 3600 self.initial_conditions[c_n] = param.n.prim.c_init_av / param.n.prim.c_max self.initial_conditions[c_p] = param.p.prim.c_init_av / param.p.prim.c_max @@ -59,9 +55,11 @@ def __init__(self, name="ECM with split OCV"): Un = param.n.prim.U(c_n, T) Up = param.p.prim.U(c_p, T) - # IR resistance, hard coded for now - IR = 0.1 - V = Up - Un - IR + # Resistance for IR expression + R = pybamm.Parameter("Ohmic resistance [Ohm]") + + # Voltage expression + V = Up - Un - I * R self.variables = { "Negative particle SOC": c_n, From dbeaeb515f690e40cbda5642f76c1fbeb78d3c96 Mon Sep 17 00:00:00 2001 From: Caitlin Parke Date: Tue, 13 Aug 2024 14:45:42 -0400 Subject: [PATCH 05/30] still working --- .../lithium_ion/basic_ecm_split_OCV.py | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py b/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py index 4cb0730958..e1e62b5b46 100644 --- a/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py +++ b/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py @@ -42,8 +42,15 @@ def __init__(self, name="ECM with split OCV"): self.initial_conditions[Q] = pybamm.Scalar(0) # Capacity in each electrode - Q_n = pybamm.Parameter("Negative electrode capacity [A.h]") - Q_p = pybamm.Parameter("Positive electrode capacity [A.h]") + Q_n = param.n.Q_init + Q_p = param.p.Q_init + # Q_n = pybamm.Parameter("Negative electrode capacity [A.h]") + # Q_p = pybamm.Parameter("Positive electrode capacity [A.h]") + R = pybamm.Parameter("Ohmic resistance [Ohm]") + Un = pybamm.Parameter("Negative electrode OCP [V]") + Up = pybamm.Parameter("Positive electrode OCP [V]") + c_n_0 = pybamm.Parameter("Negative electrode initial SOC") + c_p_0 = pybamm.Parameter("Positive electrode initial SOC") # State of charge electrode equations self.rhs[c_n] = -I / Q_n / 3600 @@ -72,6 +79,7 @@ def __init__(self, name="ECM with split OCV"): "Negative electrode potential [V]": Un, "Current variable [A]": I, "Current function [A]": I, + "Negative electrode capacity [A.h]": Q_n, } # Events specify points at which a solution should terminate From 3daabce7885a9659b71a6f79e2414d6bdb941a2f Mon Sep 17 00:00:00 2001 From: Caitlin Parke Date: Wed, 14 Aug 2024 12:01:58 -0400 Subject: [PATCH 06/30] parameter updates --- .../lithium_ion/basic_ecm_split_OCV.py | 34 ++++++++----------- 1 file changed, 15 insertions(+), 19 deletions(-) diff --git a/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py b/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py index e1e62b5b46..00fe7a2bf0 100644 --- a/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py +++ b/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py @@ -7,7 +7,7 @@ class ECMsplitOCV(BaseModel): """Basic Equivalent Circuit Model that uses two OCV functions - for each electrode from the OCV function from Lithium ion parameter sets. + for each electrode. This model is easily parameterizable with minimal parameters. This class differs from the :class: pybamm.equivalent_circuit.Thevenin() due to dual OCV functions to make up the voltage from each electrode. @@ -34,7 +34,6 @@ def __init__(self, name="ECM with split OCV"): V = pybamm.Variable("Voltage [V]") # model is isothermal - T = param.T_init I = param.current_with_time # Capacity equation @@ -42,29 +41,28 @@ def __init__(self, name="ECM with split OCV"): self.initial_conditions[Q] = pybamm.Scalar(0) # Capacity in each electrode - Q_n = param.n.Q_init - Q_p = param.p.Q_init - # Q_n = pybamm.Parameter("Negative electrode capacity [A.h]") - # Q_p = pybamm.Parameter("Positive electrode capacity [A.h]") - R = pybamm.Parameter("Ohmic resistance [Ohm]") - Un = pybamm.Parameter("Negative electrode OCP [V]") - Up = pybamm.Parameter("Positive electrode OCP [V]") - c_n_0 = pybamm.Parameter("Negative electrode initial SOC") - c_p_0 = pybamm.Parameter("Positive electrode initial SOC") + Q_n = pybamm.Parameter("Negative electrode capacity [A.h]") + Q_p = pybamm.Parameter("Positive electrode capacity [A.h]") # State of charge electrode equations + c_n_0 = pybamm.Parameter("Negative electrode initial SOC") + c_p_0 = pybamm.Parameter("Positive electrode initial SOC") self.rhs[c_n] = -I / Q_n / 3600 self.rhs[c_p] = I / Q_p / 3600 - self.initial_conditions[c_n] = param.n.prim.c_init_av / param.n.prim.c_max - self.initial_conditions[c_p] = param.p.prim.c_init_av / param.p.prim.c_max - - # OCV's for the electrodes - Un = param.n.prim.U(c_n, T) - Up = param.p.prim.U(c_p, T) + self.initial_conditions[c_n] = c_n_0 + self.initial_conditions[c_p] = c_p_0 # Resistance for IR expression R = pybamm.Parameter("Ohmic resistance [Ohm]") + # Open-circuit potential for each electrode + Un = pybamm.FunctionParameter( + "Negative electrode OCP [V]", {"Negative particle SOC": c_n} + ) + Up = pybamm.FunctionParameter( + "Positive electrode OCP [V]", {"Positive particle SOC": c_p} + ) + # Voltage expression V = Up - Un - I * R @@ -77,9 +75,7 @@ def __init__(self, name="ECM with split OCV"): "Times [s]": pybamm.t, "Positive electrode potential [V]": Up, "Negative electrode potential [V]": Un, - "Current variable [A]": I, "Current function [A]": I, - "Negative electrode capacity [A.h]": Q_n, } # Events specify points at which a solution should terminate From 758175f07874d5193467af976adf4bc67c26fd9d Mon Sep 17 00:00:00 2001 From: Caitlin Parke Date: Fri, 16 Aug 2024 00:42:04 -0400 Subject: [PATCH 07/30] added tests --- .../lithium_ion/basic_ecm_split_OCV.py | 19 +++++--- .../test_lithium_ion/test_ecm_split_OCV.py | 46 +++++++++++++++++++ .../test_lithium_ion/test_ecm_split_OCV.py | 9 ++++ 3 files changed, 67 insertions(+), 7 deletions(-) create mode 100644 tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_ecm_split_OCV.py create mode 100644 tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_ecm_split_OCV.py diff --git a/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py b/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py index 00fe7a2bf0..9195a7433d 100644 --- a/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py +++ b/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py @@ -2,10 +2,9 @@ # Equivalent Circuit Model with split OCV # import pybamm -from .base_lithium_ion_model import BaseModel -class ECMsplitOCV(BaseModel): +class ECMsplitOCV(pybamm.BaseModel): """Basic Equivalent Circuit Model that uses two OCV functions for each electrode. This model is easily parameterizable with minimal parameters. This class differs from the :class: pybamm.equivalent_circuit.Thevenin() due @@ -18,9 +17,9 @@ class ECMsplitOCV(BaseModel): """ def __init__(self, name="ECM with split OCV"): - super().__init__({}, name) + super().__init__(name) # TODO citations - param = self.param + # param = self.param ###################### # Variables @@ -34,7 +33,9 @@ def __init__(self, name="ECM with split OCV"): V = pybamm.Variable("Voltage [V]") # model is isothermal - I = param.current_with_time + I = pybamm.FunctionParameter( + "Current function [A]", {"Time [s]": pybamm.t} + ) # Capacity equation self.rhs[Q] = I / 3600 @@ -66,6 +67,10 @@ def __init__(self, name="ECM with split OCV"): # Voltage expression V = Up - Un - I * R + # Parameters for Voltage cutoff + voltage_high_cut = pybamm.Parameter('Upper voltage cut-off [V]') + voltage_low_cut = pybamm.Parameter('Lower voltage cut-off [V]') + self.variables = { "Negative particle SOC": c_n, "Positive particle SOC": c_p, @@ -80,8 +85,8 @@ def __init__(self, name="ECM with split OCV"): # Events specify points at which a solution should terminate self.events += [ - pybamm.Event("Minimum voltage [V]", V - param.voltage_low_cut), - pybamm.Event("Maximum voltage [V]", param.voltage_high_cut - V), + pybamm.Event("Minimum voltage [V]", V - voltage_low_cut), + pybamm.Event("Maximum voltage [V]", voltage_high_cut - V), pybamm.Event("Maximum Negative Electrode SOC", 0.999 - c_n), pybamm.Event("Maximum Positive Electrode SOC", 0.999 - c_p), pybamm.Event("Minimum Negative Electrode SOC", c_n - 0.0001), diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_ecm_split_OCV.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_ecm_split_OCV.py new file mode 100644 index 0000000000..e14dc7f803 --- /dev/null +++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_ecm_split_OCV.py @@ -0,0 +1,46 @@ +# +# Test that the model works should I change the base model to just regular base model from .lihtium_ion_model BaseModel +# does that set parameters that we really don't want, is this really as simple as we want it to be +import pybamm +import numpy as np + +class TestECMSplitOCVModel: + def test_run_model_with_parameters(self): + model = pybamm.lithium_ion.ECMsplitOCV() + + # example parameters + qp0 = 8.73231852 + qn0 = 5.82761507 + c0_n = 0.9013973983641687*0.9 + c0_p = 0.5142305254580026*0.83 + + # OCV functions + def Un(theta_n): + Un = 0.1493 + 0.8493*np.exp(-61.79*theta_n) + 0.3824*np.exp(-665.8*theta_n) \ + - np.exp(39.42*theta_n-41.92) - 0.03131*np.arctan(25.59*theta_n - 4.099) \ + - 0.009434*np.arctan(32.49*theta_n - 15.74) + return Un + + def Up(theta_p): + Up = -10.72*theta_p**4 + 23.88*theta_p**3 - 16.77*theta_p**2 + 2.595*theta_p + 4.563 + return Up + + pars = pybamm.ParameterValues( + { + 'Positive electrode capacity [A.h]' : qp0, + 'Ohmic resistance [Ohm]' : 0.001, + 'Negative electrode initial SOC' : c0_n, + 'Lower voltage cut-off [V]' : 2.8, + 'Positive electrode initial SOC' : c0_p, + 'Upper voltage cut-off [V]' : 4.2, + 'Negative electrode capacity [A.h]' : qn0, + 'Current function [A]' : 5, + 'Positive electrode OCP [V]' : Up, + 'Negative electrode OCP [V]' : Un, + } + ) + + # solve the model + sim = pybamm.Simulation(model, parameter_values=pars) + t_eval = np.linspace(0, 3600) + sim.solve(t_eval) \ No newline at end of file diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_ecm_split_OCV.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_ecm_split_OCV.py new file mode 100644 index 0000000000..67266de5ae --- /dev/null +++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_ecm_split_OCV.py @@ -0,0 +1,9 @@ +# +# Test for the ecm split-OCV model +# +import pybamm + +class TestECMSplitOCV: + def test_ecmsplitocv_well_posed(self): + model = pybamm.lithium_ion.ECMsplitOCV() + model.check_well_posedness() \ No newline at end of file From 99e31196650de430b775c35544e14060d87a0eb0 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 16 Aug 2024 04:42:19 +0000 Subject: [PATCH 08/30] style: pre-commit fixes --- .../lithium_ion/basic_ecm_split_OCV.py | 8 ++- .../test_lithium_ion/test_ecm_split_OCV.py | 50 ++++++++++++------- .../test_lithium_ion/test_ecm_split_OCV.py | 3 +- 3 files changed, 36 insertions(+), 25 deletions(-) diff --git a/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py b/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py index 9195a7433d..55895e6c9d 100644 --- a/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py +++ b/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py @@ -33,9 +33,7 @@ def __init__(self, name="ECM with split OCV"): V = pybamm.Variable("Voltage [V]") # model is isothermal - I = pybamm.FunctionParameter( - "Current function [A]", {"Time [s]": pybamm.t} - ) + I = pybamm.FunctionParameter("Current function [A]", {"Time [s]": pybamm.t}) # Capacity equation self.rhs[Q] = I / 3600 @@ -68,8 +66,8 @@ def __init__(self, name="ECM with split OCV"): V = Up - Un - I * R # Parameters for Voltage cutoff - voltage_high_cut = pybamm.Parameter('Upper voltage cut-off [V]') - voltage_low_cut = pybamm.Parameter('Lower voltage cut-off [V]') + voltage_high_cut = pybamm.Parameter("Upper voltage cut-off [V]") + voltage_low_cut = pybamm.Parameter("Lower voltage cut-off [V]") self.variables = { "Negative particle SOC": c_n, diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_ecm_split_OCV.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_ecm_split_OCV.py index e14dc7f803..dbf43133f3 100644 --- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_ecm_split_OCV.py +++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_ecm_split_OCV.py @@ -4,43 +4,55 @@ import pybamm import numpy as np + class TestECMSplitOCVModel: def test_run_model_with_parameters(self): model = pybamm.lithium_ion.ECMsplitOCV() - + # example parameters qp0 = 8.73231852 qn0 = 5.82761507 - c0_n = 0.9013973983641687*0.9 - c0_p = 0.5142305254580026*0.83 + c0_n = 0.9013973983641687 * 0.9 + c0_p = 0.5142305254580026 * 0.83 # OCV functions def Un(theta_n): - Un = 0.1493 + 0.8493*np.exp(-61.79*theta_n) + 0.3824*np.exp(-665.8*theta_n) \ - - np.exp(39.42*theta_n-41.92) - 0.03131*np.arctan(25.59*theta_n - 4.099) \ - - 0.009434*np.arctan(32.49*theta_n - 15.74) + Un = ( + 0.1493 + + 0.8493 * np.exp(-61.79 * theta_n) + + 0.3824 * np.exp(-665.8 * theta_n) + - np.exp(39.42 * theta_n - 41.92) + - 0.03131 * np.arctan(25.59 * theta_n - 4.099) + - 0.009434 * np.arctan(32.49 * theta_n - 15.74) + ) return Un def Up(theta_p): - Up = -10.72*theta_p**4 + 23.88*theta_p**3 - 16.77*theta_p**2 + 2.595*theta_p + 4.563 + Up = ( + -10.72 * theta_p**4 + + 23.88 * theta_p**3 + - 16.77 * theta_p**2 + + 2.595 * theta_p + + 4.563 + ) return Up - + pars = pybamm.ParameterValues( { - 'Positive electrode capacity [A.h]' : qp0, - 'Ohmic resistance [Ohm]' : 0.001, - 'Negative electrode initial SOC' : c0_n, - 'Lower voltage cut-off [V]' : 2.8, - 'Positive electrode initial SOC' : c0_p, - 'Upper voltage cut-off [V]' : 4.2, - 'Negative electrode capacity [A.h]' : qn0, - 'Current function [A]' : 5, - 'Positive electrode OCP [V]' : Up, - 'Negative electrode OCP [V]' : Un, + "Positive electrode capacity [A.h]": qp0, + "Ohmic resistance [Ohm]": 0.001, + "Negative electrode initial SOC": c0_n, + "Lower voltage cut-off [V]": 2.8, + "Positive electrode initial SOC": c0_p, + "Upper voltage cut-off [V]": 4.2, + "Negative electrode capacity [A.h]": qn0, + "Current function [A]": 5, + "Positive electrode OCP [V]": Up, + "Negative electrode OCP [V]": Un, } ) # solve the model sim = pybamm.Simulation(model, parameter_values=pars) t_eval = np.linspace(0, 3600) - sim.solve(t_eval) \ No newline at end of file + sim.solve(t_eval) diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_ecm_split_OCV.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_ecm_split_OCV.py index 67266de5ae..a26b42a459 100644 --- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_ecm_split_OCV.py +++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_ecm_split_OCV.py @@ -3,7 +3,8 @@ # import pybamm + class TestECMSplitOCV: def test_ecmsplitocv_well_posed(self): model = pybamm.lithium_ion.ECMsplitOCV() - model.check_well_posedness() \ No newline at end of file + model.check_well_posedness() From b1de9d1cf11db4e0c80817415ddca2f8d99591e9 Mon Sep 17 00:00:00 2001 From: Caitlin Parke Date: Fri, 16 Aug 2024 00:44:19 -0400 Subject: [PATCH 09/30] pre-commit --- .../lithium_ion/basic_ecm_split_OCV.py | 8 ++- .../test_lithium_ion/test_ecm_split_OCV.py | 50 ++++++++++++------- .../test_lithium_ion/test_ecm_split_OCV.py | 3 +- 3 files changed, 36 insertions(+), 25 deletions(-) diff --git a/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py b/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py index 9195a7433d..55895e6c9d 100644 --- a/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py +++ b/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py @@ -33,9 +33,7 @@ def __init__(self, name="ECM with split OCV"): V = pybamm.Variable("Voltage [V]") # model is isothermal - I = pybamm.FunctionParameter( - "Current function [A]", {"Time [s]": pybamm.t} - ) + I = pybamm.FunctionParameter("Current function [A]", {"Time [s]": pybamm.t}) # Capacity equation self.rhs[Q] = I / 3600 @@ -68,8 +66,8 @@ def __init__(self, name="ECM with split OCV"): V = Up - Un - I * R # Parameters for Voltage cutoff - voltage_high_cut = pybamm.Parameter('Upper voltage cut-off [V]') - voltage_low_cut = pybamm.Parameter('Lower voltage cut-off [V]') + voltage_high_cut = pybamm.Parameter("Upper voltage cut-off [V]") + voltage_low_cut = pybamm.Parameter("Lower voltage cut-off [V]") self.variables = { "Negative particle SOC": c_n, diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_ecm_split_OCV.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_ecm_split_OCV.py index e14dc7f803..dbf43133f3 100644 --- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_ecm_split_OCV.py +++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_ecm_split_OCV.py @@ -4,43 +4,55 @@ import pybamm import numpy as np + class TestECMSplitOCVModel: def test_run_model_with_parameters(self): model = pybamm.lithium_ion.ECMsplitOCV() - + # example parameters qp0 = 8.73231852 qn0 = 5.82761507 - c0_n = 0.9013973983641687*0.9 - c0_p = 0.5142305254580026*0.83 + c0_n = 0.9013973983641687 * 0.9 + c0_p = 0.5142305254580026 * 0.83 # OCV functions def Un(theta_n): - Un = 0.1493 + 0.8493*np.exp(-61.79*theta_n) + 0.3824*np.exp(-665.8*theta_n) \ - - np.exp(39.42*theta_n-41.92) - 0.03131*np.arctan(25.59*theta_n - 4.099) \ - - 0.009434*np.arctan(32.49*theta_n - 15.74) + Un = ( + 0.1493 + + 0.8493 * np.exp(-61.79 * theta_n) + + 0.3824 * np.exp(-665.8 * theta_n) + - np.exp(39.42 * theta_n - 41.92) + - 0.03131 * np.arctan(25.59 * theta_n - 4.099) + - 0.009434 * np.arctan(32.49 * theta_n - 15.74) + ) return Un def Up(theta_p): - Up = -10.72*theta_p**4 + 23.88*theta_p**3 - 16.77*theta_p**2 + 2.595*theta_p + 4.563 + Up = ( + -10.72 * theta_p**4 + + 23.88 * theta_p**3 + - 16.77 * theta_p**2 + + 2.595 * theta_p + + 4.563 + ) return Up - + pars = pybamm.ParameterValues( { - 'Positive electrode capacity [A.h]' : qp0, - 'Ohmic resistance [Ohm]' : 0.001, - 'Negative electrode initial SOC' : c0_n, - 'Lower voltage cut-off [V]' : 2.8, - 'Positive electrode initial SOC' : c0_p, - 'Upper voltage cut-off [V]' : 4.2, - 'Negative electrode capacity [A.h]' : qn0, - 'Current function [A]' : 5, - 'Positive electrode OCP [V]' : Up, - 'Negative electrode OCP [V]' : Un, + "Positive electrode capacity [A.h]": qp0, + "Ohmic resistance [Ohm]": 0.001, + "Negative electrode initial SOC": c0_n, + "Lower voltage cut-off [V]": 2.8, + "Positive electrode initial SOC": c0_p, + "Upper voltage cut-off [V]": 4.2, + "Negative electrode capacity [A.h]": qn0, + "Current function [A]": 5, + "Positive electrode OCP [V]": Up, + "Negative electrode OCP [V]": Un, } ) # solve the model sim = pybamm.Simulation(model, parameter_values=pars) t_eval = np.linspace(0, 3600) - sim.solve(t_eval) \ No newline at end of file + sim.solve(t_eval) diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_ecm_split_OCV.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_ecm_split_OCV.py index 67266de5ae..a26b42a459 100644 --- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_ecm_split_OCV.py +++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_ecm_split_OCV.py @@ -3,7 +3,8 @@ # import pybamm + class TestECMSplitOCV: def test_ecmsplitocv_well_posed(self): model = pybamm.lithium_ion.ECMsplitOCV() - model.check_well_posedness() \ No newline at end of file + model.check_well_posedness() From 3f1e7f50792a157ecfb0c2cb269d3f426c5032e0 Mon Sep 17 00:00:00 2001 From: Caitlin Parke Date: Mon, 19 Aug 2024 20:07:44 -0400 Subject: [PATCH 10/30] added default plotting variables --- .../lithium_ion/basic_ecm_split_OCV.py | 14 ++++++++++++-- .../test_lithium_ion/test_ecm_split_OCV.py | 4 ++-- 2 files changed, 14 insertions(+), 4 deletions(-) diff --git a/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py b/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py index 55895e6c9d..19398b4a4f 100644 --- a/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py +++ b/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py @@ -76,8 +76,8 @@ def __init__(self, name="ECM with split OCV"): "Discharge capacity [A.h]": Q, "Voltage [V]": V, "Times [s]": pybamm.t, - "Positive electrode potential [V]": Up, - "Negative electrode potential [V]": Un, + "Positive electrode OCP [V]": Up, + "Negative electrode OCP [V]": Un, "Current function [A]": I, } @@ -90,3 +90,13 @@ def __init__(self, name="ECM with split OCV"): pybamm.Event("Minimum Negative Electrode SOC", c_n - 0.0001), pybamm.Event("Minimum Positive Electrode SOC", c_p - 0.0001), ] + + @property + def default_quick_plot_variables(self): + return [ + "Voltage [V]", + ["Negative particle SOC", "Positive particle SOC"], + "Negative electrode OCP [V]", + "Positive electrode OCP [V]", + "Current [A]", + ] diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_ecm_split_OCV.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_ecm_split_OCV.py index dbf43133f3..947b3dda35 100644 --- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_ecm_split_OCV.py +++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_ecm_split_OCV.py @@ -1,6 +1,6 @@ # -# Test that the model works should I change the base model to just regular base model from .lihtium_ion_model BaseModel -# does that set parameters that we really don't want, is this really as simple as we want it to be +# Test that the model works with an example parameter set +# import pybamm import numpy as np From e4504a6d6d6a717d18326917dccb75b80d057f5e Mon Sep 17 00:00:00 2001 From: Caitlin Parke Date: Tue, 20 Aug 2024 18:18:42 -0400 Subject: [PATCH 11/30] added docs --- docs/source/api/models/lithium_ion/ecm_split_ocv.rst | 7 +++++++ docs/source/api/models/lithium_ion/index.rst | 1 + .../full_battery_models/lithium_ion/basic_ecm_split_OCV.py | 2 -- 3 files changed, 8 insertions(+), 2 deletions(-) create mode 100644 docs/source/api/models/lithium_ion/ecm_split_ocv.rst diff --git a/docs/source/api/models/lithium_ion/ecm_split_ocv.rst b/docs/source/api/models/lithium_ion/ecm_split_ocv.rst new file mode 100644 index 0000000000..2fcc4d37d8 --- /dev/null +++ b/docs/source/api/models/lithium_ion/ecm_split_ocv.rst @@ -0,0 +1,7 @@ +Equivalent Circuit Model with Split OCV (ECMsplitOCV) +===================================================== + +.. autoclass:: pybamm.lithium_ion.ECMsplitOCV + :members: + +.. footbibliography:: diff --git a/docs/source/api/models/lithium_ion/index.rst b/docs/source/api/models/lithium_ion/index.rst index 1a72c3c662..52efe44d6b 100644 --- a/docs/source/api/models/lithium_ion/index.rst +++ b/docs/source/api/models/lithium_ion/index.rst @@ -12,3 +12,4 @@ Lithium-ion Models msmr yang2017 electrode_soh + ecm_split_ocv diff --git a/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py b/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py index 19398b4a4f..2869d9b94f 100644 --- a/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py +++ b/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py @@ -18,8 +18,6 @@ class ECMsplitOCV(pybamm.BaseModel): def __init__(self, name="ECM with split OCV"): super().__init__(name) - # TODO citations - # param = self.param ###################### # Variables From 9a996c4006e2cf94ab9f39f0e96f7813f710d0bb Mon Sep 17 00:00:00 2001 From: Caitlin Parke Date: Wed, 21 Aug 2024 16:46:28 -0400 Subject: [PATCH 12/30] updated for pytest and better coverage --- .../test_lithium_ion/test_ecm_split_OCV.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_ecm_split_OCV.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_ecm_split_OCV.py index a26b42a459..5645466d93 100644 --- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_ecm_split_OCV.py +++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_ecm_split_OCV.py @@ -8,3 +8,8 @@ class TestECMSplitOCV: def test_ecmsplitocv_well_posed(self): model = pybamm.lithium_ion.ECMsplitOCV() model.check_well_posedness() + + def test_get_default_quick_plot_variables(self): + model = pybamm.lithium_ion.ECMsplitOCV() + variables = model.default_quick_plot_variables + assert "Current [A]" in variables From f23ac77d82de6ec0ea05a3da98f4873e862b75b1 Mon Sep 17 00:00:00 2001 From: Caitlin Parke Date: Tue, 3 Sep 2024 13:11:22 -0400 Subject: [PATCH 13/30] updated tests and naming convention --- CHANGELOG.md | 1 + .../api/models/lithium_ion/ecm_split_ocv.rst | 4 +-- .../lithium_ion/__init__.py | 4 +-- ...ic_ecm_split_OCV.py => basic_splitOCVR.py} | 36 +++++++++---------- ...est_ecm_split_OCV.py => test_splitOCVR.py} | 24 ++++++------- ...est_ecm_split_OCV.py => test_splitOCVR.py} | 6 ++-- 6 files changed, 37 insertions(+), 38 deletions(-) rename src/pybamm/models/full_battery_models/lithium_ion/{basic_ecm_split_OCV.py => basic_splitOCVR.py} (66%) rename tests/integration/test_models/test_full_battery_models/test_lithium_ion/{test_ecm_split_OCV.py => test_splitOCVR.py} (73%) rename tests/unit/test_models/test_full_battery_models/test_lithium_ion/{test_ecm_split_OCV.py => test_splitOCVR.py} (70%) diff --git a/CHANGELOG.md b/CHANGELOG.md index c1e8a720e7..72755d1538 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ ## Features +- Added a lithium ion equivalent circuit model with split open circuit voltages for each electrode. ([#4330](https://github.com/pybamm-team/PyBaMM/pull/4330)) - Added additional user-configurable options to the (`IDAKLUSolver`) and adjusted the default values to improve performance. ([#4282](https://github.com/pybamm-team/PyBaMM/pull/4282)) - Added the diffusion element to be used in the Thevenin model. ([#4254](https://github.com/pybamm-team/PyBaMM/pull/4254)) diff --git a/docs/source/api/models/lithium_ion/ecm_split_ocv.rst b/docs/source/api/models/lithium_ion/ecm_split_ocv.rst index 2fcc4d37d8..a7d833cf55 100644 --- a/docs/source/api/models/lithium_ion/ecm_split_ocv.rst +++ b/docs/source/api/models/lithium_ion/ecm_split_ocv.rst @@ -1,7 +1,7 @@ -Equivalent Circuit Model with Split OCV (ECMsplitOCV) +Equivalent Circuit Model with Split OCV (SplitOCVR) ===================================================== -.. autoclass:: pybamm.lithium_ion.ECMsplitOCV +.. autoclass:: pybamm.lithium_ion.SplitOCVR :members: .. footbibliography:: diff --git a/src/pybamm/models/full_battery_models/lithium_ion/__init__.py b/src/pybamm/models/full_battery_models/lithium_ion/__init__.py index f98da10b73..556e8de31c 100644 --- a/src/pybamm/models/full_battery_models/lithium_ion/__init__.py +++ b/src/pybamm/models/full_battery_models/lithium_ion/__init__.py @@ -24,9 +24,9 @@ from .Yang2017 import Yang2017 from .mpm import MPM from .msmr import MSMR -from .basic_ecm_split_OCV import ECMsplitOCV +from .basic_splitOCVR import SplitOCVR __all__ = ['Yang2017', 'base_lithium_ion_model', 'basic_dfn', 'basic_dfn_composite', 'basic_dfn_half_cell', 'basic_spm', 'dfn', 'electrode_soh', 'electrode_soh_half_cell', 'mpm', 'msmr', - 'newman_tobias', 'spm', 'spme', 'ecm_split_ocv'] + 'newman_tobias', 'spm', 'spme', 'basic_splitOCVR'] diff --git a/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py b/src/pybamm/models/full_battery_models/lithium_ion/basic_splitOCVR.py similarity index 66% rename from src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py rename to src/pybamm/models/full_battery_models/lithium_ion/basic_splitOCVR.py index 2869d9b94f..386a5c08fc 100644 --- a/src/pybamm/models/full_battery_models/lithium_ion/basic_ecm_split_OCV.py +++ b/src/pybamm/models/full_battery_models/lithium_ion/basic_splitOCVR.py @@ -4,7 +4,7 @@ import pybamm -class ECMsplitOCV(pybamm.BaseModel): +class SplitOCVR(pybamm.BaseModel): """Basic Equivalent Circuit Model that uses two OCV functions for each electrode. This model is easily parameterizable with minimal parameters. This class differs from the :class: pybamm.equivalent_circuit.Thevenin() due @@ -25,8 +25,8 @@ def __init__(self, name="ECM with split OCV"): # All variables are only time-dependent # No domain definition needed - c_n = pybamm.Variable("Negative particle SOC") - c_p = pybamm.Variable("Positive particle SOC") + theta_n = pybamm.Variable("Negative particle stoichiometry") + theta_p = pybamm.Variable("Positive particle stoichiometry") Q = pybamm.Variable("Discharge capacity [A.h]") V = pybamm.Variable("Voltage [V]") @@ -42,22 +42,22 @@ def __init__(self, name="ECM with split OCV"): Q_p = pybamm.Parameter("Positive electrode capacity [A.h]") # State of charge electrode equations - c_n_0 = pybamm.Parameter("Negative electrode initial SOC") - c_p_0 = pybamm.Parameter("Positive electrode initial SOC") - self.rhs[c_n] = -I / Q_n / 3600 - self.rhs[c_p] = I / Q_p / 3600 - self.initial_conditions[c_n] = c_n_0 - self.initial_conditions[c_p] = c_p_0 + theta_n_0 = pybamm.Parameter("Negative electrode initial stoichiometry") + theta_p_0 = pybamm.Parameter("Positive electrode initial stoichiometry") + self.rhs[theta_n] = -I / Q_n / 3600 + self.rhs[theta_p] = I / Q_p / 3600 + self.initial_conditions[theta_n] = theta_n_0 + self.initial_conditions[theta_p] = theta_p_0 # Resistance for IR expression R = pybamm.Parameter("Ohmic resistance [Ohm]") # Open-circuit potential for each electrode Un = pybamm.FunctionParameter( - "Negative electrode OCP [V]", {"Negative particle SOC": c_n} + "Negative electrode OCP [V]", {"Negative particle stoichiometry": theta_n} ) Up = pybamm.FunctionParameter( - "Positive electrode OCP [V]", {"Positive particle SOC": c_p} + "Positive electrode OCP [V]", {"Positive particle stoichiometry": theta_p} ) # Voltage expression @@ -68,8 +68,8 @@ def __init__(self, name="ECM with split OCV"): voltage_low_cut = pybamm.Parameter("Lower voltage cut-off [V]") self.variables = { - "Negative particle SOC": c_n, - "Positive particle SOC": c_p, + "Negative particle stoichiometry": theta_n, + "Positive particle stoichiometry": theta_p, "Current [A]": I, "Discharge capacity [A.h]": Q, "Voltage [V]": V, @@ -83,17 +83,17 @@ def __init__(self, name="ECM with split OCV"): self.events += [ pybamm.Event("Minimum voltage [V]", V - voltage_low_cut), pybamm.Event("Maximum voltage [V]", voltage_high_cut - V), - pybamm.Event("Maximum Negative Electrode SOC", 0.999 - c_n), - pybamm.Event("Maximum Positive Electrode SOC", 0.999 - c_p), - pybamm.Event("Minimum Negative Electrode SOC", c_n - 0.0001), - pybamm.Event("Minimum Positive Electrode SOC", c_p - 0.0001), + pybamm.Event("Maximum Negative Electrode stoichiometry", 0.999 - theta_n), + pybamm.Event("Maximum Positive Electrode stoichiometry", 0.999 - theta_p), + pybamm.Event("Minimum Negative Electrode stoichiometry", theta_n - 0.0001), + pybamm.Event("Minimum Positive Electrode stoichiometry", theta_p - 0.0001), ] @property def default_quick_plot_variables(self): return [ "Voltage [V]", - ["Negative particle SOC", "Positive particle SOC"], + ["Negative particle stoichiometry", "Positive particle stoichiometry"], "Negative electrode OCP [V]", "Positive electrode OCP [V]", "Current [A]", diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_ecm_split_OCV.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_splitOCVR.py similarity index 73% rename from tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_ecm_split_OCV.py rename to tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_splitOCVR.py index 947b3dda35..ef4391acdd 100644 --- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_ecm_split_OCV.py +++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_splitOCVR.py @@ -3,17 +3,16 @@ # import pybamm import numpy as np +import tests -class TestECMSplitOCVModel: - def test_run_model_with_parameters(self): - model = pybamm.lithium_ion.ECMsplitOCV() - +class TestSplitOCVR: + def test_basic_processing(self): # example parameters qp0 = 8.73231852 qn0 = 5.82761507 - c0_n = 0.9013973983641687 * 0.9 - c0_p = 0.5142305254580026 * 0.83 + theta0_n = 0.9013973983641687 * 0.9 + theta0_p = 0.5142305254580026 * 0.83 # OCV functions def Un(theta_n): @@ -41,18 +40,17 @@ def Up(theta_p): { "Positive electrode capacity [A.h]": qp0, "Ohmic resistance [Ohm]": 0.001, - "Negative electrode initial SOC": c0_n, + "Negative electrode initial stoichiometry": theta0_n, "Lower voltage cut-off [V]": 2.8, - "Positive electrode initial SOC": c0_p, + "Positive electrode initial stoichiometry": theta0_p, "Upper voltage cut-off [V]": 4.2, "Negative electrode capacity [A.h]": qn0, "Current function [A]": 5, "Positive electrode OCP [V]": Up, "Negative electrode OCP [V]": Un, + "Nominal cell capacity [A.h]": 5, } ) - - # solve the model - sim = pybamm.Simulation(model, parameter_values=pars) - t_eval = np.linspace(0, 3600) - sim.solve(t_eval) + model = pybamm.lithium_ion.SplitOCVR() + modeltest = tests.StandardModelTest(model) + modeltest.test_all(param=pars) diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_ecm_split_OCV.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_splitOCVR.py similarity index 70% rename from tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_ecm_split_OCV.py rename to tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_splitOCVR.py index 5645466d93..e807ec1607 100644 --- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_ecm_split_OCV.py +++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_splitOCVR.py @@ -4,12 +4,12 @@ import pybamm -class TestECMSplitOCV: +class TestSplitOCVR: def test_ecmsplitocv_well_posed(self): - model = pybamm.lithium_ion.ECMsplitOCV() + model = pybamm.lithium_ion.SplitOCVR() model.check_well_posedness() def test_get_default_quick_plot_variables(self): - model = pybamm.lithium_ion.ECMsplitOCV() + model = pybamm.lithium_ion.SplitOCVR() variables = model.default_quick_plot_variables assert "Current [A]" in variables From 3b21bf62ae604edfd322c7702550b3d97a9b9cbd Mon Sep 17 00:00:00 2001 From: "Eric G. Kratz" Date: Tue, 3 Sep 2024 16:45:19 -0400 Subject: [PATCH 14/30] Move changelog update to unreleased --- CHANGELOG.md | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 0962bdbe23..9e61a1b795 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,10 +1,13 @@ # [Unreleased](https://github.com/pybamm-team/PyBaMM/) +## Features + +- Added a lithium ion equivalent circuit model with split open circuit voltages for each electrode. ([#4330](https://github.com/pybamm-team/PyBaMM/pull/4330)) + # [v24.9.0](https://github.com/pybamm-team/PyBaMM/tree/v24.9.0) - 2024-09-03 ## Features -- Added a lithium ion equivalent circuit model with split open circuit voltages for each electrode. ([#4330](https://github.com/pybamm-team/PyBaMM/pull/4330)) - Added additional user-configurable options to the (`IDAKLUSolver`) and adjusted the default values to improve performance. ([#4282](https://github.com/pybamm-team/PyBaMM/pull/4282)) - Added the diffusion element to be used in the Thevenin model. ([#4254](https://github.com/pybamm-team/PyBaMM/pull/4254)) From 343cea4377894458a1009c6e222d7dedadf462f1 Mon Sep 17 00:00:00 2001 From: Pip Liggins Date: Fri, 13 Sep 2024 15:51:17 -0700 Subject: [PATCH 15/30] change check_extrapolation to use t/y_event --- src/pybamm/solvers/base_solver.py | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/src/pybamm/solvers/base_solver.py b/src/pybamm/solvers/base_solver.py index 1df9aef35f..fbe9ac4b04 100644 --- a/src/pybamm/solvers/base_solver.py +++ b/src/pybamm/solvers/base_solver.py @@ -864,9 +864,9 @@ def solve( # If the new initial conditions are different # and cannot be evaluated directly, set up again self.set_up(model, model_inputs_list[0], t_eval, ics_only=True) - self._model_set_up[model]["initial conditions"] = ( - model.concatenated_initial_conditions - ) + self._model_set_up[model][ + "initial conditions" + ] = model.concatenated_initial_conditions else: # Set the standard initial conditions self._set_initial_conditions(model, t_eval[0], model_inputs_list[0]) @@ -1478,8 +1478,12 @@ def check_extrapolation(self, solution, events): # second pass: check if the extrapolation events are within the tolerance last_state = solution.last_state - t = last_state.all_ts[0][0] - y = last_state.all_ys[0][:, 0] + if solution.t_event: + t = solution.t_event[0] + y = solution.y_event[:, 0] + else: + t = last_state.all_ts[0][0] + y = last_state.all_ys[0][:, 0] inputs = last_state.all_inputs[0] if isinstance(y, casadi.DM): From d6b213c41083c097a73025ba0b302140e6541ad2 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 13 Sep 2024 23:14:56 +0000 Subject: [PATCH 16/30] style: pre-commit fixes --- src/pybamm/solvers/base_solver.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/pybamm/solvers/base_solver.py b/src/pybamm/solvers/base_solver.py index fbe9ac4b04..a43a710e7b 100644 --- a/src/pybamm/solvers/base_solver.py +++ b/src/pybamm/solvers/base_solver.py @@ -864,9 +864,9 @@ def solve( # If the new initial conditions are different # and cannot be evaluated directly, set up again self.set_up(model, model_inputs_list[0], t_eval, ics_only=True) - self._model_set_up[model][ - "initial conditions" - ] = model.concatenated_initial_conditions + self._model_set_up[model]["initial conditions"] = ( + model.concatenated_initial_conditions + ) else: # Set the standard initial conditions self._set_initial_conditions(model, t_eval[0], model_inputs_list[0]) From 8e3eb31d71474b3c304d9be707ba0084329f0516 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon, 16 Sep 2024 16:00:04 -0400 Subject: [PATCH 17/30] Build(deps): bump github/codeql-action in the actions group (#4444) Bumps the actions group with 1 update: [github/codeql-action](https://github.com/github/codeql-action). Updates `github/codeql-action` from 3.26.6 to 3.26.7 - [Release notes](https://github.com/github/codeql-action/releases) - [Changelog](https://github.com/github/codeql-action/blob/main/CHANGELOG.md) - [Commits](https://github.com/github/codeql-action/compare/4dd16135b69a43b6c8efb853346f8437d92d3c93...8214744c546c1e5c8f03dde8fab3a7353211988d) --- updated-dependencies: - dependency-name: github/codeql-action dependency-type: direct:production update-type: version-update:semver-patch dependency-group: actions ... Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- .github/workflows/scorecard.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/scorecard.yml b/.github/workflows/scorecard.yml index 784ccebc4f..efc8b23141 100644 --- a/.github/workflows/scorecard.yml +++ b/.github/workflows/scorecard.yml @@ -68,6 +68,6 @@ jobs: # Upload the results to GitHub's code scanning dashboard (optional). # Commenting out will disable upload of results to your repo's Code Scanning dashboard - name: "Upload to code-scanning" - uses: github/codeql-action/upload-sarif@4dd16135b69a43b6c8efb853346f8437d92d3c93 # v3.26.6 + uses: github/codeql-action/upload-sarif@8214744c546c1e5c8f03dde8fab3a7353211988d # v3.26.7 with: sarif_file: results.sarif From 05a0b240b7208adf989ea3afe05ca9065a118570 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 16 Sep 2024 22:08:33 +0100 Subject: [PATCH 18/30] chore: update pre-commit hooks (#4445) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit updates: - [github.com/astral-sh/ruff-pre-commit: v0.6.4 → v0.6.5](https://github.com/astral-sh/ruff-pre-commit/compare/v0.6.4...v0.6.5) Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> --- .pre-commit-config.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 985cd0291a..2135fa3c71 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -4,7 +4,7 @@ ci: repos: - repo: https://github.com/astral-sh/ruff-pre-commit - rev: "v0.6.4" + rev: "v0.6.5" hooks: - id: ruff args: [--fix, --show-fixes] From e1118ecfaf18ea479d881a92fe051b7c73372d6c Mon Sep 17 00:00:00 2001 From: Marc Berliner <34451391+MarcBerliner@users.noreply.github.com> Date: Wed, 18 Sep 2024 12:35:47 -0400 Subject: [PATCH 19/30] Update CODEOWNERS (#4452) --- .github/CODEOWNERS | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/CODEOWNERS b/.github/CODEOWNERS index 0f503f09a7..6984abf32c 100644 --- a/.github/CODEOWNERS +++ b/.github/CODEOWNERS @@ -10,7 +10,7 @@ src/pybamm/meshes/ @martinjrobins @rtimms @valentinsulzer @rtimms src/pybamm/models/ @brosaplanella @DrSOKane @rtimms @valentinsulzer @TomTranter @rtimms src/pybamm/parameters/ @brosaplanella @DrSOKane @rtimms @valentinsulzer @TomTranter @rtimms @kratman src/pybamm/plotting/ @martinjrobins @rtimms @Saransh-cpp @valentinsulzer @rtimms @kratman @agriyakhetarpal -src/pybamm/solvers/ @martinjrobins @rtimms @valentinsulzer @TomTranter @rtimms +src/pybamm/solvers/ @martinjrobins @rtimms @valentinsulzer @TomTranter @rtimms @MarcBerliner src/pybamm/spatial_methods/ @martinjrobins @rtimms @valentinsulzer @rtimms src/pybamm/* @pybamm-team/maintainers # the files directly under /pybamm/, will not recurse From 48dbb68b56689fbe4e1dadd2838076fc1db68bc4 Mon Sep 17 00:00:00 2001 From: Martin Robinson Date: Wed, 18 Sep 2024 22:21:38 +0100 Subject: [PATCH 20/30] feat: add OpenMP parallelization to IDAKLU solver for lists of input parameters (#4449) * new solver option `num_solvers`, indicates how many solves run in parallel * existing `num_threads` gives total number of threads which are distributed among `num_solvers` --- CHANGELOG.md | 1 + CMakeLists.txt | 23 ++- src/pybamm/solvers/base_solver.py | 59 +++--- src/pybamm/solvers/c_solvers/idaklu.cpp | 15 +- .../solvers/c_solvers/idaklu/IDAKLUSolver.hpp | 20 +- .../c_solvers/idaklu/IDAKLUSolverGroup.cpp | 145 +++++++++++++++ .../c_solvers/idaklu/IDAKLUSolverGroup.hpp | 48 +++++ .../c_solvers/idaklu/IDAKLUSolverOpenMP.hpp | 17 +- .../c_solvers/idaklu/IDAKLUSolverOpenMP.inl | 174 ++++-------------- .../solvers/c_solvers/idaklu/Options.cpp | 17 ++ .../solvers/c_solvers/idaklu/Options.hpp | 1 + .../solvers/c_solvers/idaklu/Solution.hpp | 10 + .../solvers/c_solvers/idaklu/SolutionData.cpp | 67 +++++++ .../solvers/c_solvers/idaklu/SolutionData.hpp | 73 ++++++++ .../solvers/c_solvers/idaklu/common.cpp | 16 +- .../solvers/c_solvers/idaklu/common.hpp | 24 ++- .../c_solvers/idaklu/idaklu_solver.hpp | 122 ++++++++---- src/pybamm/solvers/idaklu_solver.py | 52 ++++-- src/pybamm/solvers/jax_solver.py | 8 +- tests/unit/test_solvers/test_idaklu_solver.py | 41 +++++ 20 files changed, 677 insertions(+), 256 deletions(-) create mode 100644 src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverGroup.cpp create mode 100644 src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverGroup.hpp create mode 100644 src/pybamm/solvers/c_solvers/idaklu/SolutionData.cpp create mode 100644 src/pybamm/solvers/c_solvers/idaklu/SolutionData.hpp diff --git a/CHANGELOG.md b/CHANGELOG.md index 80cda39db7..b65d1342f6 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ ## Features - Added sensitivity calculation support for `pybamm.Simulation` and `pybamm.Experiment` ([#4415](https://github.com/pybamm-team/PyBaMM/pull/4415)) +- Added OpenMP parallelization to IDAKLU solver for lists of input parameters ([#4449](https://github.com/pybamm-team/PyBaMM/pull/4449)) ## Optimizations - Removed the `start_step_offset` setting and disabled minimum `dt` warnings for drive cycles with the (`IDAKLUSolver`). ([#4416](https://github.com/pybamm-team/PyBaMM/pull/4416)) diff --git a/CMakeLists.txt b/CMakeLists.txt index ad56ac34ca..a7f68ce7a0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -19,7 +19,7 @@ endif() project(idaklu) -set(CMAKE_CXX_STANDARD 14) +set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_STANDARD_REQUIRED ON) set(CMAKE_CXX_EXTENSIONS OFF) set(CMAKE_EXPORT_COMPILE_COMMANDS 1) @@ -82,6 +82,8 @@ pybind11_add_module(idaklu src/pybamm/solvers/c_solvers/idaklu/idaklu_solver.hpp src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolver.cpp src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolver.hpp + src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverGroup.cpp + src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverGroup.hpp src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverOpenMP.inl src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverOpenMP.hpp src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverOpenMP_solvers.cpp @@ -94,6 +96,8 @@ pybind11_add_module(idaklu src/pybamm/solvers/c_solvers/idaklu/common.cpp src/pybamm/solvers/c_solvers/idaklu/Solution.cpp src/pybamm/solvers/c_solvers/idaklu/Solution.hpp + src/pybamm/solvers/c_solvers/idaklu/SolutionData.cpp + src/pybamm/solvers/c_solvers/idaklu/SolutionData.hpp src/pybamm/solvers/c_solvers/idaklu/Options.hpp src/pybamm/solvers/c_solvers/idaklu/Options.cpp # IDAKLU expressions / function evaluation [abstract] @@ -138,6 +142,23 @@ set_target_properties( INSTALL_RPATH_USE_LINK_PATH TRUE ) +# openmp +if (${CMAKE_SYSTEM_NAME} MATCHES "Darwin") + execute_process( + COMMAND "brew" "--prefix" + OUTPUT_VARIABLE HOMEBREW_PREFIX + OUTPUT_STRIP_TRAILING_WHITESPACE) + if (OpenMP_ROOT) + set(OpenMP_ROOT "${OpenMP_ROOT}:${HOMEBREW_PREFIX}/opt/libomp") + else() + set(OpenMP_ROOT "${HOMEBREW_PREFIX}/opt/libomp") + endif() +endif() +find_package(OpenMP) +if(OpenMP_CXX_FOUND) + target_link_libraries(idaklu PRIVATE OpenMP::OpenMP_CXX) +endif() + set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${PROJECT_SOURCE_DIR}) # Sundials find_package(SUNDIALS REQUIRED) diff --git a/src/pybamm/solvers/base_solver.py b/src/pybamm/solvers/base_solver.py index 1df9aef35f..dc1ac0cd72 100644 --- a/src/pybamm/solvers/base_solver.py +++ b/src/pybamm/solvers/base_solver.py @@ -86,6 +86,10 @@ def supports_interp(self): def root_method(self): return self._root_method + @property + def supports_parallel_solve(self): + return False + @root_method.setter def root_method(self, method): if method == "casadi": @@ -896,17 +900,8 @@ def solve( pybamm.logger.verbose( f"Calling solver for {t_eval[start_index]} < t < {t_eval[end_index - 1]}" ) - ninputs = len(model_inputs_list) - if ninputs == 1: - new_solution = self._integrate( - model, - t_eval[start_index:end_index], - model_inputs_list[0], - t_interp=t_interp, - ) - new_solutions = [new_solution] - elif model.convert_to_format == "jax": - # Jax can parallelize over the inputs efficiently + if self.supports_parallel_solve: + # Jax and IDAKLU solver can accept a list of inputs new_solutions = self._integrate( model, t_eval[start_index:end_index], @@ -914,18 +909,28 @@ def solve( t_interp, ) else: - with mp.get_context(self._mp_context).Pool(processes=nproc) as p: - new_solutions = p.starmap( - self._integrate, - zip( - [model] * ninputs, - [t_eval[start_index:end_index]] * ninputs, - model_inputs_list, - [t_interp] * ninputs, - ), + ninputs = len(model_inputs_list) + if ninputs == 1: + new_solution = self._integrate( + model, + t_eval[start_index:end_index], + model_inputs_list[0], + t_interp=t_interp, ) - p.close() - p.join() + new_solutions = [new_solution] + else: + with mp.get_context(self._mp_context).Pool(processes=nproc) as p: + new_solutions = p.starmap( + self._integrate, + zip( + [model] * ninputs, + [t_eval[start_index:end_index]] * ninputs, + model_inputs_list, + [t_interp] * ninputs, + ), + ) + p.close() + p.join() # Setting the solve time for each segment. # pybamm.Solution.__add__ assumes attribute solve_time. solve_time = timer.time() @@ -995,7 +1000,7 @@ def solve( ) # Return solution(s) - if ninputs == 1: + if len(solutions) == 1: return solutions[0] else: return solutions @@ -1350,7 +1355,13 @@ def step( # Step pybamm.logger.verbose(f"Stepping for {t_start_shifted:.0f} < t < {t_end:.0f}") timer.reset() - solution = self._integrate(model, t_eval, model_inputs, t_interp) + + # API for _integrate is different for JaxSolver and IDAKLUSolver + if self.supports_parallel_solve: + solutions = self._integrate(model, t_eval, [model_inputs], t_interp) + solution = solutions[0] + else: + solution = self._integrate(model, t_eval, model_inputs, t_interp) solution.solve_time = timer.time() # Check if extrapolation occurred diff --git a/src/pybamm/solvers/c_solvers/idaklu.cpp b/src/pybamm/solvers/c_solvers/idaklu.cpp index 3ef0194403..db7147feb2 100644 --- a/src/pybamm/solvers/c_solvers/idaklu.cpp +++ b/src/pybamm/solvers/c_solvers/idaklu.cpp @@ -9,6 +9,7 @@ #include #include "idaklu/idaklu_solver.hpp" +#include "idaklu/IDAKLUSolverGroup.hpp" #include "idaklu/IdakluJax.hpp" #include "idaklu/common.hpp" #include "idaklu/Expressions/Casadi/CasadiFunctions.hpp" @@ -26,15 +27,17 @@ casadi::Function generate_casadi_function(const std::string &data) namespace py = pybind11; PYBIND11_MAKE_OPAQUE(std::vector); +PYBIND11_MAKE_OPAQUE(std::vector); PYBIND11_MODULE(idaklu, m) { m.doc() = "sundials solvers"; // optional module docstring py::bind_vector>(m, "VectorNdArray"); + py::bind_vector>(m, "VectorSolution"); - py::class_(m, "IDAKLUSolver") - .def("solve", &IDAKLUSolver::solve, + py::class_(m, "IDAKLUSolverGroup") + .def("solve", &IDAKLUSolverGroup::solve, "perform a solve", py::arg("t_eval"), py::arg("t_interp"), @@ -43,8 +46,8 @@ PYBIND11_MODULE(idaklu, m) py::arg("inputs"), py::return_value_policy::take_ownership); - m.def("create_casadi_solver", &create_idaklu_solver, - "Create a casadi idaklu solver object", + m.def("create_casadi_solver_group", &create_idaklu_solver_group, + "Create a group of casadi idaklu solver objects", py::arg("number_of_states"), py::arg("number_of_parameters"), py::arg("rhs_alg"), @@ -70,8 +73,8 @@ PYBIND11_MODULE(idaklu, m) py::return_value_policy::take_ownership); #ifdef IREE_ENABLE - m.def("create_iree_solver", &create_idaklu_solver, - "Create a iree idaklu solver object", + m.def("create_iree_solver_group", &create_idaklu_solver_group, + "Create a group of iree idaklu solver objects", py::arg("number_of_states"), py::arg("number_of_parameters"), py::arg("rhs_alg"), diff --git a/src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolver.hpp b/src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolver.hpp index 29b451e6d3..379d64783a 100644 --- a/src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolver.hpp +++ b/src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolver.hpp @@ -2,7 +2,8 @@ #define PYBAMM_IDAKLU_CASADI_SOLVER_HPP #include "common.hpp" -#include "Solution.hpp" +#include "SolutionData.hpp" + /** * Abstract base class for solutions that can use different solvers and vector @@ -24,14 +25,17 @@ class IDAKLUSolver ~IDAKLUSolver() = default; /** - * @brief Abstract solver method that returns a Solution class + * @brief Abstract solver method that executes the solver */ - virtual Solution solve( - np_array t_eval_np, - np_array t_interp_np, - np_array y0_np, - np_array yp0_np, - np_array_dense inputs) = 0; + virtual SolutionData solve( + const std::vector &t_eval, + const std::vector &t_interp, + const realtype *y0, + const realtype *yp0, + const realtype *inputs, + bool save_adaptive_steps, + bool save_interp_steps + ) = 0; /** * Abstract method to initialize the solver, once vectors and solver classes diff --git a/src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverGroup.cpp b/src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverGroup.cpp new file mode 100644 index 0000000000..8a76d73cfe --- /dev/null +++ b/src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverGroup.cpp @@ -0,0 +1,145 @@ +#include "IDAKLUSolverGroup.hpp" +#include +#include + +std::vector IDAKLUSolverGroup::solve( + np_array t_eval_np, + np_array t_interp_np, + np_array y0_np, + np_array yp0_np, + np_array inputs) { + DEBUG("IDAKLUSolverGroup::solve"); + + // If t_interp is empty, save all adaptive steps + bool save_adaptive_steps = t_interp_np.size() == 0; + + const realtype* t_eval_begin = t_eval_np.data(); + const realtype* t_eval_end = t_eval_begin + t_eval_np.size(); + const realtype* t_interp_begin = t_interp_np.data(); + const realtype* t_interp_end = t_interp_begin + t_interp_np.size(); + + // Process the time inputs + // 1. Get the sorted and unique t_eval vector + auto const t_eval = makeSortedUnique(t_eval_begin, t_eval_end); + + // 2.1. Get the sorted and unique t_interp vector + auto const t_interp_unique_sorted = makeSortedUnique(t_interp_begin, t_interp_end); + + // 2.2 Remove the t_eval values from t_interp + auto const t_interp_setdiff = setDiff(t_interp_unique_sorted.begin(), t_interp_unique_sorted.end(), t_eval_begin, t_eval_end); + + // 2.3 Finally, get the sorted and unique t_interp vector with t_eval values removed + auto const t_interp = makeSortedUnique(t_interp_setdiff.begin(), t_interp_setdiff.end()); + + int const number_of_evals = t_eval.size(); + int const number_of_interps = t_interp.size(); + + // setDiff removes entries of t_interp that overlap with + // t_eval, so we need to check if we need to interpolate any unique points. + // This is not the same as save_adaptive_steps since some entries of t_interp + // may be removed by setDiff + bool save_interp_steps = number_of_interps > 0; + + // 3. Check if the timestepping entries are valid + if (number_of_evals < 2) { + throw std::invalid_argument( + "t_eval must have at least 2 entries" + ); + } else if (save_interp_steps) { + if (t_interp.front() < t_eval.front()) { + throw std::invalid_argument( + "t_interp values must be greater than the smallest t_eval value: " + + std::to_string(t_eval.front()) + ); + } else if (t_interp.back() > t_eval.back()) { + throw std::invalid_argument( + "t_interp values must be less than the greatest t_eval value: " + + std::to_string(t_eval.back()) + ); + } + } + + auto n_coeffs = number_of_states + number_of_parameters * number_of_states; + + // check y0 and yp0 and inputs have the correct dimensions + if (y0_np.ndim() != 2) + throw std::domain_error("y0 has wrong number of dimensions. Expected 2 but got " + std::to_string(y0_np.ndim())); + if (yp0_np.ndim() != 2) + throw std::domain_error("yp0 has wrong number of dimensions. Expected 2 but got " + std::to_string(yp0_np.ndim())); + if (inputs.ndim() != 2) + throw std::domain_error("inputs has wrong number of dimensions. Expected 2 but got " + std::to_string(inputs.ndim())); + + auto number_of_groups = y0_np.shape()[0]; + + // check y0 and yp0 and inputs have the correct shape + if (y0_np.shape()[1] != n_coeffs) + throw std::domain_error( + "y0 has wrong number of cols. Expected " + std::to_string(n_coeffs) + + " but got " + std::to_string(y0_np.shape()[1])); + + if (yp0_np.shape()[1] != n_coeffs) + throw std::domain_error( + "yp0 has wrong number of cols. Expected " + std::to_string(n_coeffs) + + " but got " + std::to_string(yp0_np.shape()[1])); + + if (yp0_np.shape()[0] != number_of_groups) + throw std::domain_error( + "yp0 has wrong number of rows. Expected " + std::to_string(number_of_groups) + + " but got " + std::to_string(yp0_np.shape()[0])); + + if (inputs.shape()[0] != number_of_groups) + throw std::domain_error( + "inputs has wrong number of rows. Expected " + std::to_string(number_of_groups) + + " but got " + std::to_string(inputs.shape()[0])); + + const std::size_t solves_per_thread = number_of_groups / m_solvers.size(); + const std::size_t remainder_solves = number_of_groups % m_solvers.size(); + + const realtype *y0 = y0_np.data(); + const realtype *yp0 = yp0_np.data(); + const realtype *inputs_data = inputs.data(); + + std::vector results(number_of_groups); + + std::optional exception; + + omp_set_num_threads(m_solvers.size()); + #pragma omp parallel for + for (int i = 0; i < m_solvers.size(); i++) { + try { + for (int j = 0; j < solves_per_thread; j++) { + const std::size_t index = i * solves_per_thread + j; + const realtype *y = y0 + index * y0_np.shape(1); + const realtype *yp = yp0 + index * yp0_np.shape(1); + const realtype *input = inputs_data + index * inputs.shape(1); + results[index] = m_solvers[i]->solve(t_eval, t_interp, y, yp, input, save_adaptive_steps, save_interp_steps); + } + } catch (std::exception &e) { + // If an exception is thrown, we need to catch it and rethrow it outside the parallel region + #pragma omp critical + { + exception = e; + } + } + } + + if (exception.has_value()) { + py::set_error(PyExc_ValueError, exception->what()); + throw py::error_already_set(); + } + + for (int i = 0; i < remainder_solves; i++) { + const std::size_t index = number_of_groups - remainder_solves + i; + const realtype *y = y0 + index * y0_np.shape(1); + const realtype *yp = yp0 + index * yp0_np.shape(1); + const realtype *input = inputs_data + index * inputs.shape(1); + results[index] = m_solvers[i]->solve(t_eval, t_interp, y, yp, input, save_adaptive_steps, save_interp_steps); + } + + // create solutions (needs to be serial as we're using the Python GIL) + std::vector solutions(number_of_groups); + for (int i = 0; i < number_of_groups; i++) { + solutions[i] = results[i].generate_solution(); + } + return solutions; +} diff --git a/src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverGroup.hpp b/src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverGroup.hpp new file mode 100644 index 0000000000..609b3b6fca --- /dev/null +++ b/src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverGroup.hpp @@ -0,0 +1,48 @@ +#ifndef PYBAMM_IDAKLU_SOLVER_GROUP_HPP +#define PYBAMM_IDAKLU_SOLVER_GROUP_HPP + +#include "IDAKLUSolver.hpp" +#include "common.hpp" + +/** + * @brief class for a group of solvers. + */ +class IDAKLUSolverGroup +{ +public: + + /** + * @brief Default constructor + */ + IDAKLUSolverGroup(std::vector> solvers, int number_of_states, int number_of_parameters): + m_solvers(std::move(solvers)), + number_of_states(number_of_states), + number_of_parameters(number_of_parameters) + {} + + // no copy constructor (unique_ptr cannot be copied) + IDAKLUSolverGroup(IDAKLUSolverGroup &) = delete; + + /** + * @brief Default destructor + */ + ~IDAKLUSolverGroup() = default; + + /** + * @brief solver method that returns a vector of Solutions + */ + std::vector solve( + np_array t_eval_np, + np_array t_interp_np, + np_array y0_np, + np_array yp0_np, + np_array inputs); + + + private: + std::vector> m_solvers; + int number_of_states; + int number_of_parameters; +}; + +#endif // PYBAMM_IDAKLU_SOLVER_GROUP_HPP diff --git a/src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverOpenMP.hpp b/src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverOpenMP.hpp index ca710fbff6..36c2872c3e 100644 --- a/src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverOpenMP.hpp +++ b/src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverOpenMP.hpp @@ -52,6 +52,7 @@ class IDAKLUSolverOpenMP : public IDAKLUSolver int const number_of_states; // cppcheck-suppress unusedStructMember int const number_of_parameters; // cppcheck-suppress unusedStructMember int const number_of_events; // cppcheck-suppress unusedStructMember + int number_of_timesteps; int precon_type; // cppcheck-suppress unusedStructMember N_Vector yy, yp, avtol; // y, y', and absolute tolerance N_Vector *yyS; // cppcheck-suppress unusedStructMember @@ -106,12 +107,16 @@ class IDAKLUSolverOpenMP : public IDAKLUSolver /** * @brief The main solve method that solves for each variable and time step */ - Solution solve( - np_array t_eval_np, - np_array t_interp_np, - np_array y0_np, - np_array yp0_np, - np_array_dense inputs) override; + SolutionData solve( + const std::vector &t_eval, + const std::vector &t_interp, + const realtype *y0, + const realtype *yp0, + const realtype *inputs, + bool save_adaptive_steps, + bool save_interp_steps + ) override; + /** * @brief Concrete implementation of initialization method diff --git a/src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverOpenMP.inl b/src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverOpenMP.inl index fd8eb38257..313c4ce12a 100644 --- a/src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverOpenMP.inl +++ b/src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverOpenMP.inl @@ -3,6 +3,7 @@ #include #include "common.hpp" +#include "SolutionData.hpp" template IDAKLUSolverOpenMP::IDAKLUSolverOpenMP( @@ -86,11 +87,20 @@ IDAKLUSolverOpenMP::IDAKLUSolverOpenMP( template void IDAKLUSolverOpenMP::AllocateVectors() { + DEBUG("IDAKLUSolverOpenMP::AllocateVectors (num_threads = " << setup_opts.num_threads << ")"); // Create vectors - yy = N_VNew_OpenMP(number_of_states, setup_opts.num_threads, sunctx); - yp = N_VNew_OpenMP(number_of_states, setup_opts.num_threads, sunctx); - avtol = N_VNew_OpenMP(number_of_states, setup_opts.num_threads, sunctx); - id = N_VNew_OpenMP(number_of_states, setup_opts.num_threads, sunctx); + if (setup_opts.num_threads == 1) { + yy = N_VNew_Serial(number_of_states, sunctx); + yp = N_VNew_Serial(number_of_states, sunctx); + avtol = N_VNew_Serial(number_of_states, sunctx); + id = N_VNew_Serial(number_of_states, sunctx); + } else { + DEBUG("IDAKLUSolverOpenMP::AllocateVectors OpenMP"); + yy = N_VNew_OpenMP(number_of_states, setup_opts.num_threads, sunctx); + yp = N_VNew_OpenMP(number_of_states, setup_opts.num_threads, sunctx); + avtol = N_VNew_OpenMP(number_of_states, setup_opts.num_threads, sunctx); + id = N_VNew_OpenMP(number_of_states, setup_opts.num_threads, sunctx); + } } template @@ -290,6 +300,7 @@ void IDAKLUSolverOpenMP::Initialize() { template IDAKLUSolverOpenMP::~IDAKLUSolverOpenMP() { + DEBUG("IDAKLUSolverOpenMP::~IDAKLUSolverOpenMP"); // Free memory if (sensitivity) { IDASensFree(ida_mem); @@ -313,63 +324,24 @@ IDAKLUSolverOpenMP::~IDAKLUSolverOpenMP() { } template -Solution IDAKLUSolverOpenMP::solve( - np_array t_eval_np, - np_array t_interp_np, - np_array y0_np, - np_array yp0_np, - np_array_dense inputs +SolutionData IDAKLUSolverOpenMP::solve( + const std::vector &t_eval, + const std::vector &t_interp, + const realtype *y0, + const realtype *yp0, + const realtype *inputs, + bool save_adaptive_steps, + bool save_interp_steps ) { DEBUG("IDAKLUSolver::solve"); + const int number_of_evals = t_eval.size(); + const int number_of_interps = t_interp.size(); - // If t_interp is empty, save all adaptive steps - bool save_adaptive_steps = t_interp_np.unchecked<1>().size() == 0; - - // Process the time inputs - // 1. Get the sorted and unique t_eval vector - auto const t_eval = makeSortedUnique(t_eval_np); - - // 2.1. Get the sorted and unique t_interp vector - auto const t_interp_unique_sorted = makeSortedUnique(t_interp_np); - - // 2.2 Remove the t_eval values from t_interp - auto const t_interp_setdiff = setDiff(t_interp_unique_sorted, t_eval); - - // 2.3 Finally, get the sorted and unique t_interp vector with t_eval values removed - auto const t_interp = makeSortedUnique(t_interp_setdiff); - - int const number_of_evals = t_eval.size(); - int const number_of_interps = t_interp.size(); - - // setDiff removes entries of t_interp that overlap with - // t_eval, so we need to check if we need to interpolate any unique points. - // This is not the same as save_adaptive_steps since some entries of t_interp - // may be removed by setDiff - bool save_interp_steps = number_of_interps > 0; - - // 3. Check if the timestepping entries are valid - if (number_of_evals < 2) { - throw std::invalid_argument( - "t_eval must have at least 2 entries" - ); - } else if (save_interp_steps) { - if (t_interp.front() < t_eval.front()) { - throw std::invalid_argument( - "t_interp values must be greater than the smallest t_eval value: " - + std::to_string(t_eval.front()) - ); - } else if (t_interp.back() > t_eval.back()) { - throw std::invalid_argument( - "t_interp values must be less than the greatest t_eval value: " - + std::to_string(t_eval.back()) - ); - } + if (t.size() < number_of_evals + number_of_interps) { + InitializeStorage(number_of_evals + number_of_interps); } - // Initialize length_of_return_vector, t, y, and yS - InitializeStorage(number_of_evals + number_of_interps); - int i_save = 0; realtype t0 = t_eval.front(); @@ -386,24 +358,11 @@ Solution IDAKLUSolverOpenMP::solve( t_interp_next = t_interp[0]; } - auto y0 = y0_np.unchecked<1>(); - auto yp0 = yp0_np.unchecked<1>(); auto n_coeffs = number_of_states + number_of_parameters * number_of_states; - if (y0.size() != n_coeffs) { - throw std::domain_error( - "y0 has wrong size. Expected " + std::to_string(n_coeffs) + - " but got " + std::to_string(y0.size())); - } else if (yp0.size() != n_coeffs) { - throw std::domain_error( - "yp0 has wrong size. Expected " + std::to_string(n_coeffs) + - " but got " + std::to_string(yp0.size())); - } - // set inputs - auto p_inputs = inputs.unchecked<2>(); for (int i = 0; i < functions->inputs.size(); i++) { - functions->inputs[i] = p_inputs(i, 0); + functions->inputs[i] = inputs[i]; } // Setup consistent initialization @@ -543,8 +502,8 @@ Solution IDAKLUSolverOpenMP::solve( PrintStats(); } - int const number_of_timesteps = i_save; - int count; + // store number of timesteps so we can generate the solution later + number_of_timesteps = i_save; // Copy the data to return as numpy arrays @@ -554,23 +513,9 @@ Solution IDAKLUSolverOpenMP::solve( t_return[i] = t[i]; } - py::capsule free_t_when_done( - t_return, - [](void *f) { - realtype *vect = reinterpret_cast(f); - delete[] vect; - } - ); - - np_array t_ret = np_array( - number_of_timesteps, - &t_return[0], - free_t_when_done - ); - // States, y realtype *y_return = new realtype[number_of_timesteps * length_of_return_vector]; - count = 0; + int count = 0; for (size_t i = 0; i < number_of_timesteps; i++) { for (size_t j = 0; j < length_of_return_vector; j++) { y_return[count] = y[i][j]; @@ -578,20 +523,6 @@ Solution IDAKLUSolverOpenMP::solve( } } - py::capsule free_y_when_done( - y_return, - [](void *f) { - realtype *vect = reinterpret_cast(f); - delete[] vect; - } - ); - - np_array y_ret = np_array( - number_of_timesteps * length_of_return_vector, - &y_return[0], - free_y_when_done - ); - // Sensitivity states, yS // Note: Ordering of vector is different if computing outputs vs returning // the complete state vector @@ -614,43 +545,7 @@ Solution IDAKLUSolverOpenMP::solve( } } - py::capsule free_yS_when_done( - yS_return, - [](void *f) { - realtype *vect = reinterpret_cast(f); - delete[] vect; - } - ); - - np_array yS_ret = np_array( - vector { - arg_sens0, - arg_sens1, - arg_sens2 - }, - &yS_return[0], - free_yS_when_done - ); - - // Final state slice, yterm - py::capsule free_yterm_when_done( - yterm_return, - [](void *f) { - realtype *vect = reinterpret_cast(f); - delete[] vect; - } - ); - - np_array y_term = np_array( - length_of_final_sv_slice, - &yterm_return[0], - free_yterm_when_done - ); - - // Store the solution - Solution sol(retval, t_ret, y_ret, yS_ret, y_term); - - return sol; + return SolutionData(retval, number_of_timesteps, length_of_return_vector, arg_sens0, arg_sens1, arg_sens2, length_of_final_sv_slice, t_return, y_return, yS_return, yterm_return); } template @@ -828,9 +723,8 @@ void IDAKLUSolverOpenMP::SetStepOutputSensitivities( template void IDAKLUSolverOpenMP::CheckErrors(int const & flag) { if (flag < 0) { - auto message = (std::string("IDA failed with flag ") + std::to_string(flag)).c_str(); - py::set_error(PyExc_ValueError, message); - throw py::error_already_set(); + auto message = std::string("IDA failed with flag ") + std::to_string(flag); + throw std::runtime_error(message.c_str()); } } diff --git a/src/pybamm/solvers/c_solvers/idaklu/Options.cpp b/src/pybamm/solvers/c_solvers/idaklu/Options.cpp index b6a33e016e..51544040ee 100644 --- a/src/pybamm/solvers/c_solvers/idaklu/Options.cpp +++ b/src/pybamm/solvers/c_solvers/idaklu/Options.cpp @@ -1,6 +1,7 @@ #include "Options.hpp" #include #include +#include using namespace std::string_literals; @@ -11,9 +12,25 @@ SetupOptions::SetupOptions(py::dict &py_opts) precon_half_bandwidth(py_opts["precon_half_bandwidth"].cast()), precon_half_bandwidth_keep(py_opts["precon_half_bandwidth_keep"].cast()), num_threads(py_opts["num_threads"].cast()), + num_solvers(py_opts["num_solvers"].cast()), linear_solver(py_opts["linear_solver"].cast()), linsol_max_iterations(py_opts["linsol_max_iterations"].cast()) { + if (num_solvers > num_threads) + { + throw std::domain_error( + "Number of solvers must be less than or equal to the number of threads" + ); + } + + // input num_threads is the overall number of threads to use. num_solvers of these + // will be used to run solvers in parallel, leaving num_threads / num_solvers threads + // to be used by each solver. From here on num_threads is the number of threads to be used by each solver + num_threads = static_cast( + std::floor( + static_cast(num_threads) / static_cast(num_solvers) + ) + ); using_sparse_matrix = true; using_banded_matrix = false; diff --git a/src/pybamm/solvers/c_solvers/idaklu/Options.hpp b/src/pybamm/solvers/c_solvers/idaklu/Options.hpp index 66a175cfff..d0c0c1d766 100644 --- a/src/pybamm/solvers/c_solvers/idaklu/Options.hpp +++ b/src/pybamm/solvers/c_solvers/idaklu/Options.hpp @@ -15,6 +15,7 @@ struct SetupOptions { int precon_half_bandwidth; int precon_half_bandwidth_keep; int num_threads; + int num_solvers; // IDALS linear solver interface std::string linear_solver; // klu, lapack, spbcg int linsol_max_iterations; diff --git a/src/pybamm/solvers/c_solvers/idaklu/Solution.hpp b/src/pybamm/solvers/c_solvers/idaklu/Solution.hpp index 72d48fa644..a43e6a7174 100644 --- a/src/pybamm/solvers/c_solvers/idaklu/Solution.hpp +++ b/src/pybamm/solvers/c_solvers/idaklu/Solution.hpp @@ -9,6 +9,11 @@ class Solution { public: + /** + * @brief Default Constructor + */ + Solution() = default; + /** * @brief Constructor */ @@ -17,6 +22,11 @@ class Solution { } + /** + * @brief Default copy from another Solution + */ + Solution(const Solution &solution) = default; + int flag; np_array t; np_array y; diff --git a/src/pybamm/solvers/c_solvers/idaklu/SolutionData.cpp b/src/pybamm/solvers/c_solvers/idaklu/SolutionData.cpp new file mode 100644 index 0000000000..00c2ddbccc --- /dev/null +++ b/src/pybamm/solvers/c_solvers/idaklu/SolutionData.cpp @@ -0,0 +1,67 @@ +#include "SolutionData.hpp" + +Solution SolutionData::generate_solution() { + py::capsule free_t_when_done( + t_return, + [](void *f) { + realtype *vect = reinterpret_cast(f); + delete[] vect; + } + ); + + np_array t_ret = np_array( + number_of_timesteps, + &t_return[0], + free_t_when_done + ); + + py::capsule free_y_when_done( + y_return, + [](void *f) { + realtype *vect = reinterpret_cast(f); + delete[] vect; + } + ); + + np_array y_ret = np_array( + number_of_timesteps * length_of_return_vector, + &y_return[0], + free_y_when_done + ); + + py::capsule free_yS_when_done( + yS_return, + [](void *f) { + realtype *vect = reinterpret_cast(f); + delete[] vect; + } + ); + + np_array yS_ret = np_array( + std::vector { + arg_sens0, + arg_sens1, + arg_sens2 + }, + &yS_return[0], + free_yS_when_done + ); + + // Final state slice, yterm + py::capsule free_yterm_when_done( + yterm_return, + [](void *f) { + realtype *vect = reinterpret_cast(f); + delete[] vect; + } + ); + + np_array y_term = np_array( + length_of_final_sv_slice, + &yterm_return[0], + free_yterm_when_done + ); + + // Store the solution + return Solution(flag, t_ret, y_ret, yS_ret, y_term); +} diff --git a/src/pybamm/solvers/c_solvers/idaklu/SolutionData.hpp b/src/pybamm/solvers/c_solvers/idaklu/SolutionData.hpp new file mode 100644 index 0000000000..815e41daca --- /dev/null +++ b/src/pybamm/solvers/c_solvers/idaklu/SolutionData.hpp @@ -0,0 +1,73 @@ +#ifndef PYBAMM_IDAKLU_SOLUTION_DATA_HPP +#define PYBAMM_IDAKLU_SOLUTION_DATA_HPP + + +#include "common.hpp" +#include "Solution.hpp" + +/** + * @brief SolutionData class. Contains all the data needed to create a Solution + */ +class SolutionData +{ + public: + /** + * @brief Default constructor + */ + SolutionData() = default; + + /** + * @brief constructor using fields + */ + SolutionData( + int flag, + int number_of_timesteps, + int length_of_return_vector, + int arg_sens0, + int arg_sens1, + int arg_sens2, + int length_of_final_sv_slice, + realtype *t_return, + realtype *y_return, + realtype *yS_return, + realtype *yterm_return): + flag(flag), + number_of_timesteps(number_of_timesteps), + length_of_return_vector(length_of_return_vector), + arg_sens0(arg_sens0), + arg_sens1(arg_sens1), + arg_sens2(arg_sens2), + length_of_final_sv_slice(length_of_final_sv_slice), + t_return(t_return), + y_return(y_return), + yS_return(yS_return), + yterm_return(yterm_return) + {} + + + /** + * @brief Default copy from another SolutionData + */ + SolutionData(const SolutionData &solution_data) = default; + + /** + * @brief Create a solution object from this data + */ + Solution generate_solution(); + +private: + + int flag; + int number_of_timesteps; + int length_of_return_vector; + int arg_sens0; + int arg_sens1; + int arg_sens2; + int length_of_final_sv_slice; + realtype *t_return; + realtype *y_return; + realtype *yS_return; + realtype *yterm_return; +}; + +#endif // PYBAMM_IDAKLU_SOLUTION_DATA_HPP diff --git a/src/pybamm/solvers/c_solvers/idaklu/common.cpp b/src/pybamm/solvers/c_solvers/idaklu/common.cpp index bf38acc56a..161c14f340 100644 --- a/src/pybamm/solvers/c_solvers/idaklu/common.cpp +++ b/src/pybamm/solvers/c_solvers/idaklu/common.cpp @@ -11,21 +11,9 @@ std::vector numpy2realtype(const np_array& input_np) { return output; } -std::vector setDiff(const std::vector& A, const std::vector& B) { - std::vector result; - if (!(A.empty())) { - std::set_difference(A.begin(), A.end(), B.begin(), B.end(), std::back_inserter(result)); - } - return result; -} -std::vector makeSortedUnique(const std::vector& input) { - std::unordered_set uniqueSet(input.begin(), input.end()); // Remove duplicates - std::vector uniqueVector(uniqueSet.begin(), uniqueSet.end()); // Convert to vector - std::sort(uniqueVector.begin(), uniqueVector.end()); // Sort the vector - return uniqueVector; -} std::vector makeSortedUnique(const np_array& input_np) { - return makeSortedUnique(numpy2realtype(input_np)); + const auto input_vec = numpy2realtype(input_np); + return makeSortedUnique(input_vec.begin(), input_vec.end()); } diff --git a/src/pybamm/solvers/c_solvers/idaklu/common.hpp b/src/pybamm/solvers/c_solvers/idaklu/common.hpp index 3289326541..58be90932e 100644 --- a/src/pybamm/solvers/c_solvers/idaklu/common.hpp +++ b/src/pybamm/solvers/c_solvers/idaklu/common.hpp @@ -31,8 +31,8 @@ #include namespace py = pybind11; -using np_array = py::array_t; -using np_array_dense = py::array_t; +// note: we rely on c_style ordering for numpy arrays so don't change this! +using np_array = py::array_t; using np_array_int = py::array_t; /** @@ -83,12 +83,25 @@ std::vector numpy2realtype(const np_array& input_np); /** * @brief Utility function to compute the set difference of two vectors */ -std::vector setDiff(const std::vector& A, const std::vector& B); +template +std::vector setDiff(const T1 a_begin, const T1 a_end, const T2 b_begin, const T2 b_end) { + std::vector result; + if (std::distance(a_begin, a_end) > 0) { + std::set_difference(a_begin, a_end, b_begin, b_end, std::back_inserter(result)); + } + return result; +} /** * @brief Utility function to make a sorted and unique vector */ -std::vector makeSortedUnique(const std::vector& input); +template +std::vector makeSortedUnique(const T input_begin, const T input_end) { + std::unordered_set uniqueSet(input_begin, input_end); // Remove duplicates + std::vector uniqueVector(uniqueSet.begin(), uniqueSet.end()); // Convert to vector + std::sort(uniqueVector.begin(), uniqueVector.end()); // Sort the vector + return uniqueVector; +} std::vector makeSortedUnique(const np_array& input_np); @@ -126,8 +139,7 @@ std::vector makeSortedUnique(const np_array& input_np); } \ std::cout << "]" << std::endl; } -#define DEBUG_v(v, M) {\ - int N = 2; \ +#define DEBUG_v(v, N) {\ std::cout << #v << "[n=" << N << "] = ["; \ for (int i = 0; i < N; i++) { \ std::cout << v[i]; \ diff --git a/src/pybamm/solvers/c_solvers/idaklu/idaklu_solver.hpp b/src/pybamm/solvers/c_solvers/idaklu/idaklu_solver.hpp index ce1765aa82..dcc1e4f8cc 100644 --- a/src/pybamm/solvers/c_solvers/idaklu/idaklu_solver.hpp +++ b/src/pybamm/solvers/c_solvers/idaklu/idaklu_solver.hpp @@ -2,6 +2,7 @@ #define PYBAMM_CREATE_IDAKLU_SOLVER_HPP #include "IDAKLUSolverOpenMP_solvers.hpp" +#include "IDAKLUSolverGroup.hpp" #include #include @@ -12,52 +13,21 @@ */ template IDAKLUSolver *create_idaklu_solver( - int number_of_states, + std::unique_ptr functions, int number_of_parameters, - const typename ExprSet::BaseFunctionType &rhs_alg, - const typename ExprSet::BaseFunctionType &jac_times_cjmass, const np_array_int &jac_times_cjmass_colptrs, const np_array_int &jac_times_cjmass_rowvals, const int jac_times_cjmass_nnz, const int jac_bandwidth_lower, const int jac_bandwidth_upper, - const typename ExprSet::BaseFunctionType &jac_action, - const typename ExprSet::BaseFunctionType &mass_action, - const typename ExprSet::BaseFunctionType &sens, - const typename ExprSet::BaseFunctionType &events, const int number_of_events, np_array rhs_alg_id, np_array atol_np, double rel_tol, int inputs_length, - const std::vector& var_fcns, - const std::vector& dvar_dy_fcns, - const std::vector& dvar_dp_fcns, - py::dict py_opts + SolverOptions solver_opts, + SetupOptions setup_opts ) { - auto setup_opts = SetupOptions(py_opts); - auto solver_opts = SolverOptions(py_opts); - auto functions = std::make_unique( - rhs_alg, - jac_times_cjmass, - jac_times_cjmass_nnz, - jac_bandwidth_lower, - jac_bandwidth_upper, - jac_times_cjmass_rowvals, - jac_times_cjmass_colptrs, - inputs_length, - jac_action, - mass_action, - sens, - events, - number_of_states, - number_of_events, - number_of_parameters, - var_fcns, - dvar_dy_fcns, - dvar_dp_fcns, - setup_opts - ); IDAKLUSolver *idakluSolver = nullptr; @@ -189,4 +159,88 @@ IDAKLUSolver *create_idaklu_solver( return idakluSolver; } +/** + * @brief Create a group of solvers using create_idaklu_solver + */ +template +IDAKLUSolverGroup *create_idaklu_solver_group( + int number_of_states, + int number_of_parameters, + const typename ExprSet::BaseFunctionType &rhs_alg, + const typename ExprSet::BaseFunctionType &jac_times_cjmass, + const np_array_int &jac_times_cjmass_colptrs, + const np_array_int &jac_times_cjmass_rowvals, + const int jac_times_cjmass_nnz, + const int jac_bandwidth_lower, + const int jac_bandwidth_upper, + const typename ExprSet::BaseFunctionType &jac_action, + const typename ExprSet::BaseFunctionType &mass_action, + const typename ExprSet::BaseFunctionType &sens, + const typename ExprSet::BaseFunctionType &events, + const int number_of_events, + np_array rhs_alg_id, + np_array atol_np, + double rel_tol, + int inputs_length, + const std::vector& var_fcns, + const std::vector& dvar_dy_fcns, + const std::vector& dvar_dp_fcns, + py::dict py_opts +) { + auto setup_opts = SetupOptions(py_opts); + auto solver_opts = SolverOptions(py_opts); + + + std::vector> solvers; + for (int i = 0; i < setup_opts.num_solvers; i++) { + // Note: we can't copy an ExprSet as it contains raw pointers to the functions + // So we create it in the loop + auto functions = std::make_unique( + rhs_alg, + jac_times_cjmass, + jac_times_cjmass_nnz, + jac_bandwidth_lower, + jac_bandwidth_upper, + jac_times_cjmass_rowvals, + jac_times_cjmass_colptrs, + inputs_length, + jac_action, + mass_action, + sens, + events, + number_of_states, + number_of_events, + number_of_parameters, + var_fcns, + dvar_dy_fcns, + dvar_dp_fcns, + setup_opts + ); + solvers.emplace_back( + std::unique_ptr( + create_idaklu_solver( + std::move(functions), + number_of_parameters, + jac_times_cjmass_colptrs, + jac_times_cjmass_rowvals, + jac_times_cjmass_nnz, + jac_bandwidth_lower, + jac_bandwidth_upper, + number_of_events, + rhs_alg_id, + atol_np, + rel_tol, + inputs_length, + solver_opts, + setup_opts + ) + ) + ); + } + + return new IDAKLUSolverGroup(std::move(solvers), number_of_states, number_of_parameters); +} + + + #endif // PYBAMM_CREATE_IDAKLU_SOLVER_HPP diff --git a/src/pybamm/solvers/idaklu_solver.py b/src/pybamm/solvers/idaklu_solver.py index 08f86b3264..ea3903b139 100644 --- a/src/pybamm/solvers/idaklu_solver.py +++ b/src/pybamm/solvers/idaklu_solver.py @@ -29,7 +29,8 @@ idaklu = importlib.util.module_from_spec(idaklu_spec) if idaklu_spec.loader: idaklu_spec.loader.exec_module(idaklu) - except ImportError: # pragma: no cover + except ImportError as e: # pragma: no cover + print(f"Error loading idaklu: {e}") idaklu_spec = None @@ -78,8 +79,10 @@ class IDAKLUSolver(pybamm.BaseSolver): options = { # Print statistics of the solver after every solve "print_stats": False, - # Number of threads available for OpenMP + # Number of threads available for OpenMP (must be greater than or equal to `num_solvers`) "num_threads": 1, + # Number of solvers to use in parallel (for solving multiple sets of input parameters in parallel) + "num_solvers": num_threads, # Evaluation engine to use for jax, can be 'jax'(native) or 'iree' "jax_evaluator": "jax", ## Linear solver interface @@ -182,6 +185,7 @@ def __init__( "precon_half_bandwidth": 5, "precon_half_bandwidth_keep": 5, "num_threads": 1, + "num_solvers": 1, "jax_evaluator": "jax", "linear_solver": "SUNLinSol_KLU", "linsol_max_iterations": 5, @@ -209,6 +213,8 @@ def __init__( if options is None: options = default_options else: + if "num_threads" in options and "num_solvers" not in options: + options["num_solvers"] = options["num_threads"] for key, value in default_options.items(): if key not in options: options[key] = value @@ -443,7 +449,7 @@ def inputs_to_dict(inputs): if model.convert_to_format == "casadi": # Serialize casadi functions - idaklu_solver_fcn = idaklu.create_casadi_solver + idaklu_solver_fcn = idaklu.create_casadi_solver_group rhs_algebraic = idaklu.generate_function(rhs_algebraic.serialize()) jac_times_cjmass = idaklu.generate_function(jac_times_cjmass.serialize()) jac_rhs_algebraic_action = idaklu.generate_function( @@ -457,7 +463,7 @@ def inputs_to_dict(inputs): and self._options["jax_evaluator"] == "iree" ): # Convert Jax functions to MLIR (also, demote to single precision) - idaklu_solver_fcn = idaklu.create_iree_solver + idaklu_solver_fcn = idaklu.create_iree_solver_group pybamm.demote_expressions_to_32bit = True if pybamm.demote_expressions_to_32bit: warnings.warn( @@ -726,7 +732,11 @@ def _check_mlir_conversion(self, name, mlir: str): def _demote_64_to_32(self, x: pybamm.EvaluatorJax): return pybamm.EvaluatorJax._demote_64_to_32(x) - def _integrate(self, model, t_eval, inputs_dict=None, t_interp=None): + @property + def supports_parallel_solve(self): + return True + + def _integrate(self, model, t_eval, inputs_list=None, t_interp=None): """ Solve a DAE model defined by residuals with initial conditions y0. @@ -736,22 +746,30 @@ def _integrate(self, model, t_eval, inputs_dict=None, t_interp=None): The model whose solution to calculate. t_eval : numeric type The times at which to stop the integration due to a discontinuity in time. - inputs_dict : dict, optional + inputs_list: list of dict, optional Any input parameters to pass to the model when solving. t_interp : None, list or ndarray, optional The times (in seconds) at which to interpolate the solution. Defaults to `None`, which returns the adaptive time-stepping times. """ - inputs_dict = inputs_dict or {} - # stack inputs - if inputs_dict: - arrays_to_stack = [np.array(x).reshape(-1, 1) for x in inputs_dict.values()] - inputs = np.vstack(arrays_to_stack) + inputs_list = inputs_list or [{}] + + # stack inputs so that they are a 2D array of shape (number_of_inputs, number_of_parameters) + if inputs_list and inputs_list[0]: + inputs = np.vstack( + [ + np.hstack([np.array(x).reshape(-1) for x in inputs_dict.values()]) + for inputs_dict in inputs_list + ] + ) else: inputs = np.array([[]]) - y0full = model.y0full - ydot0full = model.ydot0full + # stack y0full and ydot0full so they are a 2D array of shape (number_of_inputs, number_of_states + number_of_parameters * number_of_states) + # note that y0full and ydot0full are currently 1D arrays (i.e. independent of inputs), but in the future we will support + # different initial conditions for different inputs (see https://github.com/pybamm-team/PyBaMM/pull/4260). For now we just repeat the same initial conditions for each input + y0full = np.vstack([model.y0full] * len(inputs_list)) + ydot0full = np.vstack([model.ydot0full] * len(inputs_list)) atol = getattr(model, "atol", self.atol) atol = self._check_atol_type(atol, y0full.size) @@ -761,7 +779,7 @@ def _integrate(self, model, t_eval, inputs_dict=None, t_interp=None): model.convert_to_format == "jax" and self._options["jax_evaluator"] == "iree" ): - sol = self._setup["solver"].solve( + solns = self._setup["solver"].solve( t_eval, t_interp, y0full, @@ -773,6 +791,12 @@ def _integrate(self, model, t_eval, inputs_dict=None, t_interp=None): raise pybamm.SolverError("Unsupported IDAKLU solver configuration.") integration_time = timer.time() + return [ + self._post_process_solution(soln, model, integration_time, inputs_dict) + for soln, inputs_dict in zip(solns, inputs_list) + ] + + def _post_process_solution(self, sol, model, integration_time, inputs_dict): number_of_sensitivity_parameters = self._setup[ "number_of_sensitivity_parameters" ] diff --git a/src/pybamm/solvers/jax_solver.py b/src/pybamm/solvers/jax_solver.py index da5fd4983a..bfcdef1882 100644 --- a/src/pybamm/solvers/jax_solver.py +++ b/src/pybamm/solvers/jax_solver.py @@ -185,6 +185,10 @@ def solve_model_bdf(inputs): else: return jax.jit(solve_model_bdf) + @property + def supports_parallel_solve(self): + return True + def _integrate(self, model, t_eval, inputs=None, t_interp=None): """ Solve a model defined by dydt with initial conditions y0. @@ -200,7 +204,7 @@ def _integrate(self, model, t_eval, inputs=None, t_interp=None): Returns ------- - object + list of `pybamm.Solution` An object containing the times and values of the solution, as well as various diagnostic messages. @@ -301,6 +305,4 @@ async def solve_model_async(inputs_v): sol.integration_time = integration_time solutions.append(sol) - if len(solutions) == 1: - return solutions[0] return solutions diff --git a/tests/unit/test_solvers/test_idaklu_solver.py b/tests/unit/test_solvers/test_idaklu_solver.py index 32e289b3e0..213226bb4c 100644 --- a/tests/unit/test_solvers/test_idaklu_solver.py +++ b/tests/unit/test_solvers/test_idaklu_solver.py @@ -63,6 +63,47 @@ def test_ida_roberts_klu(self): true_solution = 0.1 * solution.t np.testing.assert_array_almost_equal(solution.y[0, :], true_solution) + def test_multiple_inputs(self): + model = pybamm.BaseModel() + var = pybamm.Variable("var") + rate = pybamm.InputParameter("rate") + model.rhs = {var: -rate * var} + model.initial_conditions = {var: 2} + disc = pybamm.Discretisation() + disc.process_model(model) + + for num_threads, num_solvers in [ + [1, None], + [2, None], + [8, None], + [8, 1], + [8, 2], + [8, 7], + ]: + options = {"num_threads": num_threads} + if num_solvers is not None: + options["num_solvers"] = num_solvers + solver = pybamm.IDAKLUSolver(rtol=1e-5, atol=1e-5, options=options) + t_interp = np.linspace(0, 1, 10) + t_eval = [t_interp[0], t_interp[-1]] + ninputs = 8 + inputs_list = [{"rate": 0.01 * (i + 1)} for i in range(ninputs)] + + solutions = solver.solve( + model, t_eval, inputs=inputs_list, t_interp=t_interp + ) + + # check solution + for inputs, solution in zip(inputs_list, solutions): + print("checking solution", inputs, solution.all_inputs) + np.testing.assert_array_equal(solution.t, t_interp) + np.testing.assert_allclose( + solution.y[0], + 2 * np.exp(-inputs["rate"] * solution.t), + atol=1e-4, + rtol=1e-4, + ) + def test_model_events(self): for form in ["casadi", "iree"]: if (form == "iree") and (not pybamm.has_jax() or not pybamm.has_iree()): From 37c94f9972e7dd889593589f1a51820468e843db Mon Sep 17 00:00:00 2001 From: Pip Liggins Date: Thu, 19 Sep 2024 10:10:32 -0700 Subject: [PATCH 21/30] Add test --- tests/unit/test_solvers/test_idaklu_solver.py | 20 +++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/tests/unit/test_solvers/test_idaklu_solver.py b/tests/unit/test_solvers/test_idaklu_solver.py index 213226bb4c..045c7b6aec 100644 --- a/tests/unit/test_solvers/test_idaklu_solver.py +++ b/tests/unit/test_solvers/test_idaklu_solver.py @@ -1162,3 +1162,23 @@ def test_python_idaklu_deprecation_errors(self): match="Unsupported evaluation engine for convert_to_format=jax", ): _ = solver.solve(model, t_eval) + + def test_extrapolation_events_with_output_variables(self): + # Make sure the extrapolation checks work with output variables + model = pybamm.BaseModel() + v = pybamm.Variable("v") + c = pybamm.Variable("c") + model.variables = {"v": v, "c": c} + model.rhs = {v: -1, c: 0} + model.initial_conditions = {v: 1, c: 2} + model.events.append( + pybamm.Event( + "Ignored event", + v + 10, + pybamm.EventType.INTERPOLANT_EXTRAPOLATION, + ) + ) + solver = pybamm.IDAKLUSolver(output_variables=["c"]) + solver.set_up(model) + + solver.solve(model, t_eval=[0, 1]) From cee2cad836b28f35de0bc7fb2e9fae17f3c6e1d8 Mon Sep 17 00:00:00 2001 From: Pip Liggins Date: Thu, 19 Sep 2024 10:22:21 -0700 Subject: [PATCH 22/30] update changelog --- CHANGELOG.md | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ebe4abdd66..18bf754946 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,15 +4,16 @@ - Added sensitivity calculation support for `pybamm.Simulation` and `pybamm.Experiment` ([#4415](https://github.com/pybamm-team/PyBaMM/pull/4415)) - Added OpenMP parallelization to IDAKLU solver for lists of input parameters ([#4449](https://github.com/pybamm-team/PyBaMM/pull/4449)) +- Added phase-dependent particle options to LAM #4369 - Added a lithium ion equivalent circuit model with split open circuit voltages for each electrode (`SplitOCVR`). ([#4330](https://github.com/pybamm-team/PyBaMM/pull/4330)) ## Optimizations - Removed the `start_step_offset` setting and disabled minimum `dt` warnings for drive cycles with the (`IDAKLUSolver`). ([#4416](https://github.com/pybamm-team/PyBaMM/pull/4416)) -## Features +## Bug Fixes -- Added phase-dependent particle options to LAM #4369 +- Fixed bug where IDAKLU solver failed when `output variables` were specified and an extrapolation event is present. ([#4440](https://github.com/pybamm-team/PyBaMM/pull/4440)) ## Breaking changes From 7be637cd409508e990f04939414ed0cdea6732c2 Mon Sep 17 00:00:00 2001 From: Marc Berliner <34451391+MarcBerliner@users.noreply.github.com> Date: Thu, 19 Sep 2024 15:39:27 -0400 Subject: [PATCH 23/30] Faster (re)initialization of ODEs in `IDA` (#4453) * Fast ODE reinitialization * Update CHANGELOG.md * remove redundant `IDAReInit` call * address comments * Update IDAKLUSolverOpenMP.inl --- CHANGELOG.md | 1 + .../c_solvers/idaklu/IDAKLUSolverOpenMP.hpp | 29 ++++++- .../c_solvers/idaklu/IDAKLUSolverOpenMP.inl | 81 +++++++++++++++---- src/pybamm/solvers/idaklu_solver.py | 2 +- 4 files changed, 95 insertions(+), 18 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ebe4abdd66..9ad756a9e5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,7 @@ ## Optimizations +- Improved performance of initialization and reinitialization of ODEs in the (`IDAKLUSolver`). ([#4453](https://github.com/pybamm-team/PyBaMM/pull/4453)) - Removed the `start_step_offset` setting and disabled minimum `dt` warnings for drive cycles with the (`IDAKLUSolver`). ([#4416](https://github.com/pybamm-team/PyBaMM/pull/4416)) ## Features diff --git a/src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverOpenMP.hpp b/src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverOpenMP.hpp index 36c2872c3e..92eede3643 100644 --- a/src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverOpenMP.hpp +++ b/src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverOpenMP.hpp @@ -54,7 +54,7 @@ class IDAKLUSolverOpenMP : public IDAKLUSolver int const number_of_events; // cppcheck-suppress unusedStructMember int number_of_timesteps; int precon_type; // cppcheck-suppress unusedStructMember - N_Vector yy, yp, avtol; // y, y', and absolute tolerance + N_Vector yy, yp, y_cache, avtol; // y, y', y cache vector, and absolute tolerance N_Vector *yyS; // cppcheck-suppress unusedStructMember N_Vector *ypS; // cppcheck-suppress unusedStructMember N_Vector id; // rhs_alg_id @@ -70,6 +70,7 @@ class IDAKLUSolverOpenMP : public IDAKLUSolver vector res_dvar_dp; bool const sensitivity; // cppcheck-suppress unusedStructMember bool const save_outputs_only; // cppcheck-suppress unusedStructMember + bool is_ODE; // cppcheck-suppress unusedStructMember int length_of_return_vector; // cppcheck-suppress unusedStructMember vector t; // cppcheck-suppress unusedStructMember vector> y; // cppcheck-suppress unusedStructMember @@ -158,6 +159,32 @@ class IDAKLUSolverOpenMP : public IDAKLUSolver */ void PrintStats(); + /** + * @brief Set a consistent initialization for ODEs + */ + void ReinitializeIntegrator(const realtype& t_val); + + /** + * @brief Set a consistent initialization for the system of equations + */ + void ConsistentInitialization( + const realtype& t_val, + const realtype& t_next, + const int& icopt); + + /** + * @brief Set a consistent initialization for DAEs + */ + void ConsistentInitializationDAE( + const realtype& t_val, + const realtype& t_next, + const int& icopt); + + /** + * @brief Set a consistent initialization for ODEs + */ + void ConsistentInitializationODE(const realtype& t_val); + /** * @brief Extend the adaptive arrays by 1 */ diff --git a/src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverOpenMP.inl b/src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverOpenMP.inl index 313c4ce12a..56f546facf 100644 --- a/src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverOpenMP.inl +++ b/src/pybamm/solvers/c_solvers/idaklu/IDAKLUSolverOpenMP.inl @@ -83,6 +83,10 @@ IDAKLUSolverOpenMP::IDAKLUSolverOpenMP( if (this->setup_opts.preconditioner != "none") { precon_type = SUN_PREC_LEFT; } + + // The default is to solve a DAE for generality. This may be changed + // to an ODE during the Initialize() call + is_ODE = false; } template @@ -92,12 +96,14 @@ void IDAKLUSolverOpenMP::AllocateVectors() { if (setup_opts.num_threads == 1) { yy = N_VNew_Serial(number_of_states, sunctx); yp = N_VNew_Serial(number_of_states, sunctx); + y_cache = N_VNew_Serial(number_of_states, sunctx); avtol = N_VNew_Serial(number_of_states, sunctx); id = N_VNew_Serial(number_of_states, sunctx); } else { DEBUG("IDAKLUSolverOpenMP::AllocateVectors OpenMP"); yy = N_VNew_OpenMP(number_of_states, setup_opts.num_threads, sunctx); yp = N_VNew_OpenMP(number_of_states, setup_opts.num_threads, sunctx); + y_cache = N_VNew_OpenMP(number_of_states, setup_opts.num_threads, sunctx); avtol = N_VNew_OpenMP(number_of_states, setup_opts.num_threads, sunctx); id = N_VNew_OpenMP(number_of_states, setup_opts.num_threads, sunctx); } @@ -289,9 +295,13 @@ void IDAKLUSolverOpenMP::Initialize() { realtype *id_val; id_val = N_VGetArrayPointer(id); - int ii; - for (ii = 0; ii < number_of_states; ii++) { + // Determine if the system is an ODE + is_ODE = number_of_states > 0; + for (int ii = 0; ii < number_of_states; ii++) { id_val[ii] = id_np_val[ii]; + // check if id_val[ii] approximately equals 1 (>0.999) handles + // cases where id_val[ii] is not exactly 1 due to numerical errors + is_ODE &= id_val[ii] > 0.999; } // Variable types: differential (1) and algebraic (0) @@ -312,6 +322,7 @@ IDAKLUSolverOpenMP::~IDAKLUSolverOpenMP() { N_VDestroy(avtol); N_VDestroy(yy); N_VDestroy(yp); + N_VDestroy(y_cache); N_VDestroy(id); if (sensitivity) { @@ -386,21 +397,16 @@ SolutionData IDAKLUSolverOpenMP::solve( SetSolverOptions(); - CheckErrors(IDAReInit(ida_mem, t0, yy, yp)); - if (sensitivity) { - CheckErrors(IDASensReInit(ida_mem, IDA_SIMULTANEOUS, yyS, ypS)); - } - // Prepare first time step i_eval = 1; realtype t_eval_next = t_eval[i_eval]; + // Consistent initialization + ReinitializeIntegrator(t0); int const init_type = solver_opts.init_all_y_ic ? IDA_Y_INIT : IDA_YA_YDP_INIT; if (solver_opts.calc_ic) { - DEBUG("IDACalcIC"); - // IDACalcIC will throw a warning if it fails to find initial conditions - IDACalcIC(ida_mem, init_type, t_eval_next); + ConsistentInitialization(t0, t_eval_next, init_type); } if (sensitivity) { @@ -480,12 +486,8 @@ SolutionData IDAKLUSolverOpenMP::solve( CheckErrors(IDASetStopTime(ida_mem, t_eval_next)); // Reinitialize the solver to deal with the discontinuity at t = t_val. - // We must reinitialize the algebraic terms, so do not use init_type. - IDACalcIC(ida_mem, IDA_YA_YDP_INIT, t_eval_next); - CheckErrors(IDAReInit(ida_mem, t_val, yy, yp)); - if (sensitivity) { - CheckErrors(IDASensReInit(ida_mem, IDA_SIMULTANEOUS, yyS, ypS)); - } + ReinitializeIntegrator(t_val); + ConsistentInitialization(t_val, t_eval_next, IDA_YA_YDP_INIT); } t_prev = t_val; @@ -563,6 +565,53 @@ void IDAKLUSolverOpenMP::ExtendAdaptiveArrays() { } } +template +void IDAKLUSolverOpenMP::ReinitializeIntegrator(const realtype& t_val) { + DEBUG("IDAKLUSolver::ReinitializeIntegrator"); + CheckErrors(IDAReInit(ida_mem, t_val, yy, yp)); + if (sensitivity) { + CheckErrors(IDASensReInit(ida_mem, IDA_SIMULTANEOUS, yyS, ypS)); + } +} + +template +void IDAKLUSolverOpenMP::ConsistentInitialization( + const realtype& t_val, + const realtype& t_next, + const int& icopt) { + DEBUG("IDAKLUSolver::ConsistentInitialization"); + + if (is_ODE && icopt == IDA_YA_YDP_INIT) { + ConsistentInitializationODE(t_val); + } else { + ConsistentInitializationDAE(t_val, t_next, icopt); + } +} + +template +void IDAKLUSolverOpenMP::ConsistentInitializationDAE( + const realtype& t_val, + const realtype& t_next, + const int& icopt) { + DEBUG("IDAKLUSolver::ConsistentInitializationDAE"); + IDACalcIC(ida_mem, icopt, t_next); +} + +template +void IDAKLUSolverOpenMP::ConsistentInitializationODE( + const realtype& t_val) { + DEBUG("IDAKLUSolver::ConsistentInitializationODE"); + + // For ODEs where the mass matrix M = I, we can simplify the problem + // by analytically computing the yp values. If we take our implicit + // DAE system res(t,y,yp) = f(t,y) - I*yp, then yp = res(t,y,0). This + // avoids an expensive call to IDACalcIC. + realtype *y_cache_val = N_VGetArrayPointer(y_cache); + std::memset(y_cache_val, 0, number_of_states * sizeof(realtype)); + // Overwrite yp + residual_eval(t_val, yy, y_cache, yp, functions.get()); +} + template void IDAKLUSolverOpenMP::SetStep( realtype &tval, diff --git a/src/pybamm/solvers/idaklu_solver.py b/src/pybamm/solvers/idaklu_solver.py index ea3903b139..32048d89c0 100644 --- a/src/pybamm/solvers/idaklu_solver.py +++ b/src/pybamm/solvers/idaklu_solver.py @@ -988,7 +988,7 @@ def _rhs_dot_consistent_initialization(self, y0, model, time, inputs_dict): rhs0 = rhs_alg0[: model.len_rhs] - # for the differential terms, ydot = -M^-1 * (rhs) + # for the differential terms, ydot = M^-1 * (rhs) ydot0[: model.len_rhs] = model.mass_matrix_inv.entries @ rhs0 return ydot0 From b01268462a7b0453817819137c01e9acfa5fcf9e Mon Sep 17 00:00:00 2001 From: Pip Liggins Date: Thu, 19 Sep 2024 17:16:35 -0700 Subject: [PATCH 24/30] Switch test to triggered event --- tests/unit/test_solvers/test_idaklu_solver.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/tests/unit/test_solvers/test_idaklu_solver.py b/tests/unit/test_solvers/test_idaklu_solver.py index 045c7b6aec..b049729ae3 100644 --- a/tests/unit/test_solvers/test_idaklu_solver.py +++ b/tests/unit/test_solvers/test_idaklu_solver.py @@ -1173,12 +1173,13 @@ def test_extrapolation_events_with_output_variables(self): model.initial_conditions = {v: 1, c: 2} model.events.append( pybamm.Event( - "Ignored event", - v + 10, + "Triggered event", + v - 0.5, pybamm.EventType.INTERPOLANT_EXTRAPOLATION, ) ) solver = pybamm.IDAKLUSolver(output_variables=["c"]) solver.set_up(model) - solver.solve(model, t_eval=[0, 1]) + with pytest.warns(pybamm.SolverWarning, match="extrapolation occurred for"): + solver.solve(model, t_eval=[0, 1]) From 18ef2b20e9c9cc7e692076b4274911bc3f39433b Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon, 23 Sep 2024 16:05:52 -0400 Subject: [PATCH 25/30] Build(deps): bump github/codeql-action in the actions group (#4460) Bumps the actions group with 1 update: [github/codeql-action](https://github.com/github/codeql-action). Updates `github/codeql-action` from 3.26.7 to 3.26.8 - [Release notes](https://github.com/github/codeql-action/releases) - [Changelog](https://github.com/github/codeql-action/blob/main/CHANGELOG.md) - [Commits](https://github.com/github/codeql-action/compare/8214744c546c1e5c8f03dde8fab3a7353211988d...294a9d92911152fe08befb9ec03e240add280cb3) --- updated-dependencies: - dependency-name: github/codeql-action dependency-type: direct:production update-type: version-update:semver-patch dependency-group: actions ... Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- .github/workflows/scorecard.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/scorecard.yml b/.github/workflows/scorecard.yml index efc8b23141..3186f0e49b 100644 --- a/.github/workflows/scorecard.yml +++ b/.github/workflows/scorecard.yml @@ -68,6 +68,6 @@ jobs: # Upload the results to GitHub's code scanning dashboard (optional). # Commenting out will disable upload of results to your repo's Code Scanning dashboard - name: "Upload to code-scanning" - uses: github/codeql-action/upload-sarif@8214744c546c1e5c8f03dde8fab3a7353211988d # v3.26.7 + uses: github/codeql-action/upload-sarif@294a9d92911152fe08befb9ec03e240add280cb3 # v3.26.8 with: sarif_file: results.sarif From 1b6ef033646d7d89a0bebe12732a1d84ffee6739 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 23 Sep 2024 17:19:44 -0400 Subject: [PATCH 26/30] chore: update pre-commit hooks (#4461) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit updates: - [github.com/astral-sh/ruff-pre-commit: v0.6.5 → v0.6.7](https://github.com/astral-sh/ruff-pre-commit/compare/v0.6.5...v0.6.7) Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> --- .pre-commit-config.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 2135fa3c71..306118e254 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -4,7 +4,7 @@ ci: repos: - repo: https://github.com/astral-sh/ruff-pre-commit - rev: "v0.6.5" + rev: "v0.6.7" hooks: - id: ruff args: [--fix, --show-fixes] From e148006fc46f75edbddeec07ee943434793daa62 Mon Sep 17 00:00:00 2001 From: Abhishek Chaudhari <91185083+AbhishekChaudharii@users.noreply.github.com> Date: Wed, 25 Sep 2024 22:01:34 +0530 Subject: [PATCH 27/30] Added __iter__ and __contains__ method to pybamm.ParameterValues (#4465) * Added __contains__ and __iter__ in pybamm.ParameterValues * After running pre_commit * added tests to verify __iter__ and __contains__ method work-2 * Update tests/unit/test_parameters/test_parameter_values.py Deleted print statement which was intended for debugging Co-authored-by: Agriya Khetarpal <74401230+agriyakhetarpal@users.noreply.github.com> * Update tests/unit/test_parameters/test_parameter_values.py --------- Co-authored-by: Agriya Khetarpal <74401230+agriyakhetarpal@users.noreply.github.com> Co-authored-by: Arjun Verma --- src/pybamm/parameters/parameter_values.py | 6 ++++++ .../test_parameters/test_parameter_values.py | 21 +++++++++++++++++++ 2 files changed, 27 insertions(+) diff --git a/src/pybamm/parameters/parameter_values.py b/src/pybamm/parameters/parameter_values.py index 43c1ea17ce..d5d5878cf2 100644 --- a/src/pybamm/parameters/parameter_values.py +++ b/src/pybamm/parameters/parameter_values.py @@ -930,3 +930,9 @@ def print_evaluated_parameters(self, evaluated_parameters, output_file): file.write((s + " : {:10.4g}\n").format(name, value)) else: file.write((s + " : {:10.3E}\n").format(name, value)) + + def __contains__(self, key): + return key in self._dict_items + + def __iter__(self): + return iter(self._dict_items) diff --git a/tests/unit/test_parameters/test_parameter_values.py b/tests/unit/test_parameters/test_parameter_values.py index 4086abea6d..543d3a6b4c 100644 --- a/tests/unit/test_parameters/test_parameter_values.py +++ b/tests/unit/test_parameters/test_parameter_values.py @@ -17,6 +17,7 @@ ) from pybamm.expression_tree.exceptions import OptionError import casadi +from pybamm.parameters.parameter_values import ParameterValues class TestParameterValues: @@ -1029,3 +1030,23 @@ def test_exchange_current_density_plating(self): match="referring to the reaction at the surface of a lithium metal electrode", ): parameter_values.evaluate(param) + + def test_contains_method(self): + """Test for __contains__ method to check the functionality of 'in' keyword""" + parameter_values = ParameterValues( + {"Negative particle radius [m]": 1e-6, "Positive particle radius [m]": 2e-6} + ) + assert ( + "Negative particle radius [m]" in parameter_values + ), "Key should be found in parameter_values" + assert ( + "Invalid key" not in parameter_values + ), "Non-existent key should not be found" + + def test_iter_method(self): + """Test for __iter__ method to check if we can iterate over keys""" + parameter_values = ParameterValues( + values={"Negative particle radius [m]": 1e-6} + ) + pv = [i for i in parameter_values] + assert len(pv) == 5, "Should have 5 keys" From ffad85c766cb3c12f051e5688f469fea17e4cf2e Mon Sep 17 00:00:00 2001 From: "allcontributors[bot]" <46447321+allcontributors[bot]@users.noreply.github.com> Date: Wed, 25 Sep 2024 17:38:22 +0100 Subject: [PATCH 28/30] docs: add AbhishekChaudharii as a contributor for test (#4467) * docs: update all_contributors.md [skip ci] * docs: update README.md [skip ci] * docs: update .all-contributorsrc [skip ci] --------- Co-authored-by: allcontributors[bot] <46447321+allcontributors[bot]@users.noreply.github.com> Co-authored-by: Eric G. Kratz --- .all-contributorsrc | 3 ++- all_contributors.md | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/.all-contributorsrc b/.all-contributorsrc index 00b86f3474..caadc74b78 100644 --- a/.all-contributorsrc +++ b/.all-contributorsrc @@ -747,7 +747,8 @@ "profile": "https://github.com/AbhishekChaudharii", "contributions": [ "doc", - "code" + "code", + "test" ] }, { diff --git a/all_contributors.md b/all_contributors.md index 49a5f1d1ef..b701ef2a13 100644 --- a/all_contributors.md +++ b/all_contributors.md @@ -91,7 +91,7 @@ Thanks goes to these wonderful people ([emoji key](https://allcontributors.org/d Agnik Bakshi
Agnik Bakshi
📖 RuiheLi
RuiheLi
💻 ⚠️ chmabaur
chmabaur
🐛 💻 - Abhishek Chaudhari
Abhishek Chaudhari
📖 💻 + Abhishek Chaudhari
Abhishek Chaudhari
📖 💻 ⚠️ Shubham Bhardwaj
Shubham Bhardwaj
🚇 Jonathan Lauber
Jonathan Lauber
🚇 From 1390ea306e4200e04e6d63a9eac42333a75961e8 Mon Sep 17 00:00:00 2001 From: "Eric G. Kratz" Date: Wed, 25 Sep 2024 13:31:21 -0400 Subject: [PATCH 29/30] Remove chemistry dep warning (#4466) * Remove chemistry dep warning * Fix changelog * style: pre-commit fixes --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> --- CHANGELOG.md | 5 ++++- src/pybamm/parameters/parameter_values.py | 10 +--------- tests/unit/test_parameters/test_parameter_values.py | 6 ------ 3 files changed, 5 insertions(+), 16 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index cef03a025e..7070ebaf52 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,7 +4,8 @@ - Added sensitivity calculation support for `pybamm.Simulation` and `pybamm.Experiment` ([#4415](https://github.com/pybamm-team/PyBaMM/pull/4415)) - Added OpenMP parallelization to IDAKLU solver for lists of input parameters ([#4449](https://github.com/pybamm-team/PyBaMM/pull/4449)) -- Added phase-dependent particle options to LAM #4369 +- Added phase-dependent particle options to LAM + ([#4369](https://github.com/pybamm-team/PyBaMM/pull/4369)) - Added a lithium ion equivalent circuit model with split open circuit voltages for each electrode (`SplitOCVR`). ([#4330](https://github.com/pybamm-team/PyBaMM/pull/4330)) ## Optimizations @@ -18,6 +19,8 @@ ## Breaking changes +- Removed the deprecation warning for the chemistry argument in + ParameterValues ([#4466](https://github.com/pybamm-team/PyBaMM/pull/4466)) - The parameters "... electrode OCP entropic change [V.K-1]" and "... electrode volume change" are now expected to be functions of stoichiometry only instead of functions of both stoichiometry and maximum concentration ([#4427](https://github.com/pybamm-team/PyBaMM/pull/4427)) - Renamed `set_events` function to `add_events_from` to better reflect its purpose. ([#4421](https://github.com/pybamm-team/PyBaMM/pull/4421)) diff --git a/src/pybamm/parameters/parameter_values.py b/src/pybamm/parameters/parameter_values.py index d5d5878cf2..0a0e49cd8f 100644 --- a/src/pybamm/parameters/parameter_values.py +++ b/src/pybamm/parameters/parameter_values.py @@ -35,15 +35,7 @@ class ParameterValues: """ - def __init__(self, values, chemistry=None): - if chemistry is not None: - raise ValueError( - "The 'chemistry' keyword argument has been deprecated. " - "Call `ParameterValues` with a dictionary dictionary of " - "parameter values, or the name of a parameter set (string), " - "as the single argument, e.g. `ParameterValues('Chen2020')`.", - ) - + def __init__(self, values): # add physical constants as default values self._dict_items = pybamm.FuzzyDict( { diff --git a/tests/unit/test_parameters/test_parameter_values.py b/tests/unit/test_parameters/test_parameter_values.py index 543d3a6b4c..cdcfb30ede 100644 --- a/tests/unit/test_parameters/test_parameter_values.py +++ b/tests/unit/test_parameters/test_parameter_values.py @@ -37,12 +37,6 @@ def test_init(self): param = pybamm.ParameterValues({"a": 1, "chemistry": "lithium-ion"}) assert "chemistry" not in param.keys() - # chemistry kwarg removed - with pytest.raises( - ValueError, match="'chemistry' keyword argument has been deprecated" - ): - pybamm.ParameterValues(None, chemistry="lithium-ion") - # junk param values rejected with pytest.raises(ValueError, match="'Junk' is not a valid parameter set."): pybamm.ParameterValues("Junk") From 62a7ee8bdc2aee2e71a80fd05876b28dc541310d Mon Sep 17 00:00:00 2001 From: Brady Planden <55357039+BradyPlanden@users.noreply.github.com> Date: Thu, 26 Sep 2024 11:48:26 +0100 Subject: [PATCH 30/30] perf: refactor and speed-ups for Jax BDF Solver (#4456) * Performance refactor for Jax BDF, "BDF" as default for JaxSolver, bugfixes for calculate_sensitivities, adds JAX vectorised example * update docstring, add changelog entry * feat: adds property for explicit sensitivity attribute, suggested changes from review * examples: adds JIT compiled comparison * Apply suggestions from code review Co-authored-by: Martin Robinson * fix: post suggestions property alignment * tests: adds calculate_sensitivities check for JaxSolver * tests: move calculate_senstivities unit test for coverage --------- Co-authored-by: Martin Robinson --- CHANGELOG.md | 1 + examples/scripts/multiprocess_jax_solver.py | 57 +++ noxfile.py | 2 +- src/pybamm/solvers/base_solver.py | 12 +- src/pybamm/solvers/idaklu_solver.py | 4 + src/pybamm/solvers/jax_bdf_solver.py | 397 +++++++++----------- src/pybamm/solvers/jax_solver.py | 12 +- tests/unit/test_solvers/test_base_solver.py | 1 + tests/unit/test_solvers/test_jax_solver.py | 8 + 9 files changed, 262 insertions(+), 232 deletions(-) create mode 100644 examples/scripts/multiprocess_jax_solver.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 7070ebaf52..f70aae0272 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,6 +10,7 @@ ## Optimizations +- Performance refactor of JAX BDF Solver with default Jax method set to `"BDF"`. ([#4456](https://github.com/pybamm-team/PyBaMM/pull/4456)) - Improved performance of initialization and reinitialization of ODEs in the (`IDAKLUSolver`). ([#4453](https://github.com/pybamm-team/PyBaMM/pull/4453)) - Removed the `start_step_offset` setting and disabled minimum `dt` warnings for drive cycles with the (`IDAKLUSolver`). ([#4416](https://github.com/pybamm-team/PyBaMM/pull/4416)) diff --git a/examples/scripts/multiprocess_jax_solver.py b/examples/scripts/multiprocess_jax_solver.py new file mode 100644 index 0000000000..8192256ed1 --- /dev/null +++ b/examples/scripts/multiprocess_jax_solver.py @@ -0,0 +1,57 @@ +import pybamm +import time +import numpy as np + + +# This script provides an example for massively vectorised +# model solves using the JAX BDF solver. First, +# we set up the model and process parameters +model = pybamm.lithium_ion.SPM() +model.convert_to_format = "jax" +model.events = [] # remove events (not supported in jax) +geometry = model.default_geometry +param = pybamm.ParameterValues("Chen2020") +param.update({"Current function [A]": "[input]"}) +param.process_geometry(geometry) +param.process_model(model) + +# Discretise and setup solver +mesh = pybamm.Mesh(geometry, model.default_submesh_types, model.default_var_pts) +disc = pybamm.Discretisation(mesh, model.default_spatial_methods) +disc.process_model(model) +t_eval = np.linspace(0, 3600, 100) +solver = pybamm.JaxSolver(atol=1e-6, rtol=1e-6, method="BDF") + +# Set number of vectorised solves +values = np.linspace(0.01, 1.0, 1000) +inputs = [{"Current function [A]": value} for value in values] + +# Run solve for all inputs, with a just-in-time compilation +# occurring on the first solve. All sequential solves will +# use the compiled code, with a large performance improvement. +start_time = time.time() +sol = solver.solve(model, t_eval, inputs=inputs) +print(f"Time taken: {time.time() - start_time}") # 1.3s + +# Rerun the vectorised solve, showing performance improvement +start_time = time.time() +compiled_sol = solver.solve(model, t_eval, inputs=inputs) +print(f"Compiled time taken: {time.time() - start_time}") # 0.42s + +# Plot one of the solves +plot = pybamm.QuickPlot( + sol[5], + [ + "Negative particle concentration [mol.m-3]", + "Electrolyte concentration [mol.m-3]", + "Positive particle concentration [mol.m-3]", + "Current [A]", + "Negative electrode potential [V]", + "Electrolyte potential [V]", + "Positive electrode potential [V]", + "Voltage [V]", + ], + time_unit="seconds", + spatial_unit="um", +) +plot.dynamic_plot() diff --git a/noxfile.py b/noxfile.py index 6567ed167c..14bafcca47 100644 --- a/noxfile.py +++ b/noxfile.py @@ -222,7 +222,7 @@ def run_scripts(session): # https://bitbucket.org/pybtex-devs/pybtex/issues/169/replace-pkg_resources-with # is fixed session.install("setuptools", silent=False) - session.install("-e", ".[all,dev]", silent=False) + session.install("-e", ".[all,dev,jax]", silent=False) session.run("python", "-m", "pytest", "-m", "scripts") diff --git a/src/pybamm/solvers/base_solver.py b/src/pybamm/solvers/base_solver.py index 6334169cf0..a7fc79fe3a 100644 --- a/src/pybamm/solvers/base_solver.py +++ b/src/pybamm/solvers/base_solver.py @@ -90,6 +90,10 @@ def root_method(self): def supports_parallel_solve(self): return False + @property + def requires_explicit_sensitivities(self): + return True + @root_method.setter def root_method(self, method): if method == "casadi": @@ -141,7 +145,7 @@ def set_up(self, model, inputs=None, t_eval=None, ics_only=False): # see if we need to form the explicit sensitivity equations calculate_sensitivities_explicit = ( - model.calculate_sensitivities and not isinstance(self, pybamm.IDAKLUSolver) + model.calculate_sensitivities and self.requires_explicit_sensitivities ) self._set_up_model_sensitivities_inplace( @@ -494,11 +498,7 @@ def _set_up_model_sensitivities_inplace( # if we have a mass matrix, we need to extend it def extend_mass_matrix(M): M_extend = [M.entries] * (num_parameters + 1) - M_extend_pybamm = pybamm.Matrix(block_diag(M_extend, format="csr")) - return M_extend_pybamm - - model.mass_matrix = extend_mass_matrix(model.mass_matrix) - model.mass_matrix = extend_mass_matrix(model.mass_matrix) + return pybamm.Matrix(block_diag(M_extend, format="csr")) model.mass_matrix = extend_mass_matrix(model.mass_matrix) diff --git a/src/pybamm/solvers/idaklu_solver.py b/src/pybamm/solvers/idaklu_solver.py index 32048d89c0..2b2d852697 100644 --- a/src/pybamm/solvers/idaklu_solver.py +++ b/src/pybamm/solvers/idaklu_solver.py @@ -736,6 +736,10 @@ def _demote_64_to_32(self, x: pybamm.EvaluatorJax): def supports_parallel_solve(self): return True + @property + def requires_explicit_sensitivities(self): + return False + def _integrate(self, model, t_eval, inputs_list=None, t_interp=None): """ Solve a DAE model defined by residuals with initial conditions y0. diff --git a/src/pybamm/solvers/jax_bdf_solver.py b/src/pybamm/solvers/jax_bdf_solver.py index 6f0c62b9a8..3db82ca0da 100644 --- a/src/pybamm/solvers/jax_bdf_solver.py +++ b/src/pybamm/solvers/jax_bdf_solver.py @@ -111,68 +111,63 @@ def fun_bind_inputs(y, t): jac_bind_inputs = jax.jacfwd(fun_bind_inputs, argnums=0) - t0 = t_eval[0] - h0 = t_eval[1] - t0 - + t0, h0 = t_eval[0], t_eval[1] - t_eval[0] stepper = _bdf_init( fun_bind_inputs, jac_bind_inputs, mass, t0, y0, h0, rtol, atol ) - i = 0 y_out = jnp.empty((len(t_eval), len(y0)), dtype=y0.dtype) - init_state = [stepper, t_eval, i, y_out] - def cond_fun(state): - _, t_eval, i, _ = state + _, _, i, _ = state return i < len(t_eval) def body_fun(state): stepper, t_eval, i, y_out = state stepper = _bdf_step(stepper, fun_bind_inputs, jac_bind_inputs) - index = jnp.searchsorted(t_eval, stepper.t) - index = index.astype( - "int" + t_eval.dtype.name[-2:] - ) # Coerce index to correct type + index = jnp.searchsorted(t_eval, stepper.t).astype(jnp.int32) - def for_body(j, y_out): - t = t_eval[j] - y_out = y_out.at[jnp.index_exp[j, :]].set(_bdf_interpolate(stepper, t)) - return y_out + def interpolate_and_update(j, y_out): + y = _bdf_interpolate(stepper, t_eval[j]) + return y_out.at[j].set(y) - y_out = jax.lax.fori_loop(i, index, for_body, y_out) - return [stepper, t_eval, index, y_out] + y_out = jax.lax.fori_loop(i, index, interpolate_and_update, y_out) + return stepper, t_eval, index, y_out + + init_state = (stepper, t_eval, 0, y_out) + _, _, _, y_out = jax.lax.while_loop(cond_fun, body_fun, init_state) - stepper, t_eval, i, y_out = jax.lax.while_loop(cond_fun, body_fun, init_state) return y_out - BDFInternalStates = [ - "t", - "atol", - "rtol", - "M", - "newton_tol", - "order", - "h", - "n_equal_steps", - "D", - "y0", - "scale_y0", - "kappa", - "gamma", - "alpha", - "c", - "error_const", - "J", - "LU", - "U", - "psi", - "n_function_evals", - "n_jacobian_evals", - "n_lu_decompositions", - "n_steps", - "consistent_y0_failed", - ] - BDFState = collections.namedtuple("BDFState", BDFInternalStates) + BDFState = collections.namedtuple( + "BDFState", + [ + "t", + "atol", + "rtol", + "M", + "newton_tol", + "order", + "h", + "n_equal_steps", + "D", + "y0", + "scale_y0", + "kappa", + "gamma", + "alpha", + "c", + "error_const", + "J", + "LU", + "U", + "psi", + "n_function_evals", + "n_jacobian_evals", + "n_lu_decompositions", + "n_steps", + "consistent_y0_failed", + ], + ) jax.tree_util.register_pytree_node( BDFState, lambda xs: (tuple(xs), None), lambda _, xs: BDFState(*xs) @@ -211,62 +206,70 @@ def _bdf_init(fun, jac, mass, t0, y0, h0, rtol, atol): absolute tolerance for the solver """ - state = {} - state["t"] = t0 - state["atol"] = atol - state["rtol"] = rtol - state["M"] = mass EPS = jnp.finfo(y0.dtype).eps - state["newton_tol"] = jnp.maximum(10 * EPS / rtol, jnp.minimum(0.03, rtol**0.5)) + # Scaling and tolerance initialisation scale_y0 = atol + rtol * jnp.abs(y0) + newton_tol = jnp.maximum(10 * EPS / rtol, jnp.minimum(0.03, rtol**0.5)) + y0, not_converged = _select_initial_conditions( - fun, mass, t0, y0, state["newton_tol"], scale_y0 + fun, mass, t0, y0, newton_tol, scale_y0 ) - state["consistent_y0_failed"] = not_converged + # Compute initial function and step size f0 = fun(y0, t0) - order = 1 - state["order"] = order - state["h"] = _select_initial_step(atol, rtol, fun, t0, y0, f0, h0) - state["n_equal_steps"] = 0 + h = _select_initial_step(atol, rtol, fun, t0, y0, f0, h0) + + # Initialise the difference matrix, D D = jnp.empty((MAX_ORDER + 1, len(y0)), dtype=y0.dtype) D = D.at[jnp.index_exp[0, :]].set(y0) - D = D.at[jnp.index_exp[1, :]].set(f0 * state["h"]) - state["D"] = D - state["y0"] = y0 - state["scale_y0"] = scale_y0 + D = D.at[jnp.index_exp[1, :]].set(f0 * h) # kappa values for difference orders, taken from Table 1 of [1] kappa = jnp.array([0, -0.1850, -1 / 9, -0.0823, -0.0415, 0]) gamma = jnp.hstack((0, jnp.cumsum(1 / jnp.arange(1, MAX_ORDER + 1)))) alpha = 1.0 / ((1 - kappa) * gamma) - c = state["h"] * alpha[order] + c = h * alpha[1] error_const = kappa * gamma + 1 / jnp.arange(1, MAX_ORDER + 2) - state["kappa"] = kappa - state["gamma"] = gamma - state["alpha"] = alpha - state["c"] = c - state["error_const"] = error_const - + # Jacobian and LU decomp J = jac(y0, t0) - state["J"] = J - - state["LU"] = jax.scipy.linalg.lu_factor(state["M"] - c * J) - - state["U"] = _compute_R(order, 1) - state["psi"] = None - - state["n_function_evals"] = 2 - state["n_jacobian_evals"] = 1 - state["n_lu_decompositions"] = 1 - state["n_steps"] = 0 + LU = jax.scipy.linalg.lu_factor(mass - c * J) + U = _compute_R(1, 1) # Order 1 + + # Create initial BDFState + state = BDFState( + t=t0, + atol=atol, + rtol=rtol, + M=mass, + newton_tol=newton_tol, + consistent_y0_failed=not_converged, + order=1, + h=h, + n_equal_steps=0, + D=D, + y0=y0, + scale_y0=scale_y0, + kappa=kappa, + gamma=gamma, + alpha=alpha, + c=c, + error_const=error_const, + J=J, + LU=LU, + U=U, + psi=None, + n_function_evals=2, + n_jacobian_evals=1, + n_lu_decompositions=1, + n_steps=0, + ) - tuple_state = BDFState(*[state[k] for k in BDFInternalStates]) - y0, scale_y0 = _predict(tuple_state, D) - psi = _update_psi(tuple_state, D) - return tuple_state._replace(y0=y0, scale_y0=scale_y0, psi=psi) + # Predict initial y0, scale_yo, update state + y0, scale_y0 = _predict(state, D) + psi = _update_psi(state, D) + return state._replace(y0=y0, scale_y0=scale_y0, psi=psi) def _compute_R(order, factor): """ @@ -374,10 +377,8 @@ def _predict(state, D): """ predict forward to new step (eq 2 in [1]) """ - n = len(state.y0) - order = state.order - orders = jnp.repeat(jnp.arange(MAX_ORDER + 1).reshape(-1, 1), n, axis=1) - subD = jnp.where(orders <= order, D, 0) + orders = jnp.arange(MAX_ORDER + 1)[:, None] + subD = jnp.where(orders <= state.order, D, 0) y0 = jnp.sum(subD, axis=0) scale_y0 = state.atol + state.rtol * jnp.abs(state.y0) return y0, scale_y0 @@ -397,7 +398,7 @@ def _update_psi(state, D): def _update_difference_for_next_step(state, d): """ - update of difference equations can be done efficiently + Update of difference equations can be done efficiently by reusing d and D. From first equation on page 4 of [1]: @@ -409,34 +410,21 @@ def _update_difference_for_next_step(state, d): Combining these gives the following algorithm """ order = state.order - D = state.D - D = D.at[jnp.index_exp[order + 2]].set(d - D[order + 1]) - D = D.at[jnp.index_exp[order + 1]].set(d) - i = order - while_state = [i, D] - - def while_cond(while_state): - i, _ = while_state - return i >= 0 + D = state.D.at[order + 2].set(d - state.D[order + 1]) + D = D.at[order + 1].set(d) - def while_body(while_state): - i, D = while_state - D = D.at[jnp.index_exp[i]].add(D[i + 1]) - i -= 1 - return [i, D] + def update_D(i, D): + return D.at[order - i].add(D[order - i + 1]) - i, D = jax.lax.while_loop(while_cond, while_body, while_state) - - return D + return jax.lax.fori_loop(0, order + 1, update_D, D) def _update_step_size_and_lu(state, factor): + """ + Update step size and recompute LU decomposition. + """ state = _update_step_size(state, factor) - - # redo lu (c has changed) LU = jax.scipy.linalg.lu_factor(state.M - state.c * state.J) - n_lu_decompositions = state.n_lu_decompositions + 1 - - return state._replace(LU=LU, n_lu_decompositions=n_lu_decompositions) + return state._replace(LU=LU, n_lu_decompositions=state.n_lu_decompositions + 1) def _update_step_size(state, factor): """ @@ -449,7 +437,6 @@ def _update_step_size(state, factor): """ order = state.order h = state.h * factor - n_equal_steps = 0 c = h * state.alpha[order] # update D using equations in section 3.2 of [1] @@ -461,19 +448,14 @@ def _update_step_size(state, factor): RU = jnp.where( jnp.logical_and(I <= order, J <= order), RU, jnp.identity(MAX_ORDER + 1) ) - D = state.D - D = jnp.dot(RU.T, D) - # D = jax.ops.index_update(D, jax.ops.index[:order + 1], - # jnp.dot(RU.T, D[:order + 1])) + D = jnp.dot(RU.T, state.D) - # update psi (D has changed) + # update psi, y0 (D has changed) psi = _update_psi(state, D) - - # update y0 (D has changed) y0, scale_y0 = _predict(state, D) return state._replace( - n_equal_steps=n_equal_steps, + n_equal_steps=0, h=h, c=c, D=D, @@ -484,27 +466,23 @@ def _update_step_size(state, factor): def _update_jacobian(state, jac): """ - we update the jacobian using J(t_{n+1}, y^0_{n+1}) + Update the jacobian using J(t_{n+1}, y^0_{n+1}) following the scipy bdf implementation rather than J(t_n, y_n) as per [1] """ J = jac(state.y0, state.t + state.h) - n_jacobian_evals = state.n_jacobian_evals + 1 LU = jax.scipy.linalg.lu_factor(state.M - state.c * J) - n_lu_decompositions = state.n_lu_decompositions + 1 return state._replace( J=J, - n_jacobian_evals=n_jacobian_evals, + n_jacobian_evals=state.n_jacobian_evals + 1, LU=LU, - n_lu_decompositions=n_lu_decompositions, + n_lu_decompositions=state.n_lu_decompositions + 1, ) def _newton_iteration(state, fun): - tol = state.newton_tol - c = state.c - psi = state.psi + """ + Perform Newton iteration to solve the system. + """ y0 = state.y0 - LU = state.LU - M = state.M scale_y0 = state.scale_y0 t = state.t + state.h d = jnp.zeros(y0.shape, dtype=y0.dtype) @@ -522,17 +500,20 @@ def while_cond(while_state): def while_body(while_state): k, converged, dy_norm_old, d, y, n_function_evals = while_state + f_eval = fun(y, t) n_function_evals += 1 - b = c * f_eval - M @ (psi + d) - dy = jax.scipy.linalg.lu_solve(LU, b) + b = state.c * f_eval - state.M @ (state.psi + d) + dy = jax.scipy.linalg.lu_solve(state.LU, b) dy_norm = jnp.sqrt(jnp.mean((dy / scale_y0) ** 2)) rate = dy_norm / dy_norm_old # if iteration is not going to converge in NEWTON_MAXITER # (assuming the current rate), then abort pred = rate >= 1 - pred += rate ** (NEWTON_MAXITER - k) / (1 - rate) * dy_norm > tol + pred += ( + rate ** (NEWTON_MAXITER - k) / (1 - rate) * dy_norm > state.newton_tol + ) pred *= dy_norm_old >= 0 k += pred * (NEWTON_MAXITER - k - 1) @@ -541,7 +522,7 @@ def while_body(while_state): # if converged then break out of iteration early pred = dy_norm_old >= 0.0 - pred *= rate / (1 - rate) * dy_norm < tol + pred *= rate / (1 - rate) * dy_norm < state.newton_tol converged = (dy_norm == 0.0) + pred dy_norm_old = dy_norm @@ -564,7 +545,6 @@ def _prepare_next_step(state, d): def _prepare_next_step_order_change(state, d, y, n_iter): order = state.order - D = _update_difference_for_next_step(state, d) # Note: we are recalculating these from the while loop above, could re-use? @@ -586,7 +566,6 @@ def _prepare_next_step_order_change(state, d, y, n_iter): rms_norm(state.error_const[order + 1] * D[order + 2] / scale_y), jnp.inf, ) - error_norms = jnp.array([error_m_norm, error_norm, error_p_norm]) factors = error_norms ** (-1 / (jnp.arange(3) + order)) @@ -595,111 +574,89 @@ def _prepare_next_step_order_change(state, d, y, n_iter): max_index = jnp.argmax(factors) order += max_index - 1 + # New step size factor factor = jnp.minimum(MAX_FACTOR, safety * factors[max_index]) - - new_state = _update_step_size_and_lu(state._replace(D=D, order=order), factor) - return new_state + new_state = state._replace(D=D, order=order) + return _update_step_size_and_lu(new_state, factor) def _bdf_step(state, fun, jac): - # print('bdf_step', state.t, state.h) - # we will try and use the old jacobian unless convergence of newton iteration - # fails - updated_jacobian = False - # initialise step size and try to make the step, - # iterate, reducing step size until error is in bounds - step_accepted = False - y = jnp.empty_like(state.y0) - d = jnp.empty_like(state.y0) - n_iter = -1 - - # loop until step is accepted - while_state = [state, step_accepted, updated_jacobian, y, d, n_iter] + """ + Perform a BDF step. - def while_cond(while_state): - _, step_accepted, _, _, _, _ = while_state - return step_accepted == False # noqa: E712 + We will try and use the old jacobian unless + convergence of newton iteration fails. + """ - def while_body(while_state): - state, step_accepted, updated_jacobian, y, d, n_iter = while_state + def step_iteration(while_state): + state, updated_jacobian = while_state - # solve BDF equation using y0 as starting point + # Solve BDF equation using Newton iteration converged, n_iter, y, d, state = _newton_iteration(state, fun) - not_converged = converged == False # noqa: E712 - - # newton iteration did not converge, but jacobian has already been - # evaluated so reduce step size by 0.3 (as per [1]) and try again - state = tree_map( - partial(jnp.where, not_converged * updated_jacobian), - _update_step_size_and_lu(state, 0.3), - state, - ) - # if not_converged * updated_jacobian: - # print('not converged, update step size by 0.3') - # if not_converged * (updated_jacobian == False): - # print('not converged, update jacobian') - - # if not converged and jacobian not updated, then update the jacobian and - # try again - (state, updated_jacobian) = tree_map( - partial( - jnp.where, - not_converged * (updated_jacobian == False), # noqa: E712 + # Update Jacobian or reduce step size if not converged + # Evaluated so reduce step size by 0.3 (as per [1]) and try again + state, updated_jacobian = jax.lax.cond( + ~converged, + lambda s, uj: jax.lax.cond( + uj, + lambda s: (_update_step_size_and_lu(s, 0.3), True), + lambda s: (_update_jacobian(s, jac), True), + s, ), - (_update_jacobian(state, jac), True), - (state, False + updated_jacobian), + lambda s, uj: (s, uj), + state, + updated_jacobian, ) safety = 0.9 * (2 * NEWTON_MAXITER + 1) / (2 * NEWTON_MAXITER + n_iter) scale_y = state.atol + state.rtol * jnp.abs(y) + # Calculate error and updated step size factor # combine eq 3, 4 and 6 from [1] to obtain error # Note that error = C_k * h^{k+1} y^{k+1} # and d = D^{k+1} y_{n+1} \approx h^{k+1} y^{k+1} error = state.error_const[state.order] * d - error_norm = rms_norm(error / scale_y) - # calculate optimal step size factor as per eq 2.46 of [2] - factor = jnp.maximum( - MIN_FACTOR, safety * error_norm ** (-1 / (state.order + 1)) + # Calculate optimal step size factor as per eq 2.46 of [2] + factor = jnp.clip( + safety * error_norm ** (-1 / (state.order + 1)), MIN_FACTOR, None ) - # if converged * (error_norm > 1): - # print( - # "converged, but error is too large", - # error_norm, - # factor, - # d, - # scale_y, - # ) - - (state, step_accepted) = tree_map( - partial(jnp.where, converged * (error_norm > 1)), - (_update_step_size_and_lu(state, factor), False), - (state, converged), + # Update step size if error is too large + state = jax.lax.cond( + converged & (error_norm > 1), + lambda s: _update_step_size_and_lu(s, factor), + lambda s: s, + state, ) - return [state, step_accepted, updated_jacobian, y, d, n_iter] - - state, step_accepted, updated_jacobian, y, d, n_iter = jax.lax.while_loop( - while_cond, while_body, while_state + step_accepted = converged & (error_norm <= 1) + return (state, updated_jacobian), (step_accepted, y, d, n_iter) + + # Iterate until step is accepted + (state, _), (_, y, d, n_iter) = jax.lax.while_loop( + lambda carry_and_aux: ~carry_and_aux[1][0], + lambda carry_and_aux: step_iteration(carry_and_aux[0]), + ( + (state, False), + (False, jnp.empty_like(state.y0), jnp.empty_like(state.y0), -1), + ), ) - # take the accepted step + # Update state for the accepted step n_steps = state.n_steps + 1 t = state.t + state.h - - # a change in order is only done after running at order k for k + 1 steps - # (see page 83 of [2]) n_equal_steps = state.n_equal_steps + 1 - state = state._replace(n_equal_steps=n_equal_steps, t=t, n_steps=n_steps) - state = tree_map( - partial(jnp.where, n_equal_steps < state.order + 1), - _prepare_next_step(state, d), - _prepare_next_step_order_change(state, d, y, n_iter), + # Prepare for the next step, potentially changing order + # (see page 83 of [2]) + state = jax.lax.cond( + n_equal_steps < state.order + 1, + lambda s: _prepare_next_step(s, d), + lambda s: _prepare_next_step_order_change(s, d, y, n_iter), + state, ) return state @@ -710,8 +667,6 @@ def _bdf_interpolate(state, t_eval): definition of the interpolating polynomial can be found on page 7 of [1] """ - order = state.order - t = state.t h = state.h D = state.D j = 0 @@ -721,11 +676,11 @@ def _bdf_interpolate(state, t_eval): def while_cond(while_state): j, _, _ = while_state - return j < order + return j < state.order def while_body(while_state): j, time_factor, order_summation = while_state - time_factor *= (t_eval - (t - h * j)) / (h * (1 + j)) + time_factor *= (t_eval - (state.t - h * j)) / (h * (1 + j)) order_summation += D[j + 1] * time_factor j += 1 return [j, time_factor, order_summation] @@ -972,13 +927,13 @@ def jax_bdf_integrate(func, y0, t_eval, *args, rtol=1e-6, atol=1e-6, mass=None): """ Backward Difference formula (BDF) implicit multistep integrator. The basic algorithm is derived in :footcite:t:`byrne1975polyalgorithm`. This particular implementation - follows that implemented in the Matlab routine ode15s described in - :footcite:t:`shampine1997matlab` and the SciPy implementation - :footcite:t:`Virtanen2020` which features the NDF formulas for improved stability, - with associated differences in the error constants, and calculates the jacobian at - J(t_{n+1}, y^0_{n+1}). This implementation was based on that implemented in the - SciPy library :footcite:t:`Virtanen2020`, which also mainly follows - :footcite:t:`shampine1997matlab` but uses the more standard jacobian update. + follows the Matlab routine ode15s described in :footcite:t:`shampine1997matlab` + and the SciPy implementation :footcite:t:`Virtanen2020` which features + the NDF formulas for improved stability, with associated differences in the + error constants, and calculates the jacobian at J(t_{n+1}, y^0_{n+1}). This + implementation was based on that implemented in the SciPy library + :footcite:t:`Virtanen2020`, which also mainly follows :footcite:t:`shampine1997matlab` + but uses the more standard jacobian update. Parameters ---------- diff --git a/src/pybamm/solvers/jax_solver.py b/src/pybamm/solvers/jax_solver.py index bfcdef1882..a1f1733ed6 100644 --- a/src/pybamm/solvers/jax_solver.py +++ b/src/pybamm/solvers/jax_solver.py @@ -31,10 +31,10 @@ class JaxSolver(pybamm.BaseSolver): Parameters ---------- method: str, optional (see `jax.experimental.ode.odeint` for details) - * 'RK45' (default) uses jax.experimental.ode.odeint - * 'BDF' uses custom jax_bdf_integrate (see `jax_bdf_integrate.py` for details) + * 'BDF' (default) uses custom jax_bdf_integrate (see `jax_bdf_integrate.py` for details) + * 'RK45' uses jax.experimental.ode.odeint root_method: str, optional - Method to use to calculate consistent initial conditions. By default this uses + Method to use to calculate consistent initial conditions. By default, this uses the newton chord method internal to the jax bdf solver, otherwise choose from the set of default options defined in docs for pybamm.BaseSolver rtol : float, optional @@ -52,7 +52,7 @@ class JaxSolver(pybamm.BaseSolver): def __init__( self, - method="RK45", + method="BDF", root_method=None, rtol=1e-6, atol=1e-6, @@ -189,6 +189,10 @@ def solve_model_bdf(inputs): def supports_parallel_solve(self): return True + @property + def requires_explicit_sensitivities(self): + return False + def _integrate(self, model, t_eval, inputs=None, t_interp=None): """ Solve a model defined by dydt with initial conditions y0. diff --git a/tests/unit/test_solvers/test_base_solver.py b/tests/unit/test_solvers/test_base_solver.py index 6753513e72..1733cafc4c 100644 --- a/tests/unit/test_solvers/test_base_solver.py +++ b/tests/unit/test_solvers/test_base_solver.py @@ -19,6 +19,7 @@ def test_base_solver_init(self): assert solver.rtol == 1e-5 solver.rtol = 1e-7 assert solver.rtol == 1e-7 + assert solver.requires_explicit_sensitivities def test_root_method_init(self): solver = pybamm.BaseSolver(root_method="casadi") diff --git a/tests/unit/test_solvers/test_jax_solver.py b/tests/unit/test_solvers/test_jax_solver.py index f7e5b8d3b6..8f43eda3c7 100644 --- a/tests/unit/test_solvers/test_jax_solver.py +++ b/tests/unit/test_solvers/test_jax_solver.py @@ -237,3 +237,11 @@ def test_get_solve(self): y = solver({"rate": 0.2}) np.testing.assert_allclose(y[0], np.exp(-0.2 * t_eval), rtol=1e-6, atol=1e-6) + + # Reset solver, test passing `calculate_sensitivities` + for method in ["RK45", "BDF"]: + solver = pybamm.JaxSolver(method=method, rtol=1e-8, atol=1e-8) + solution_sens = solver.solve( + model, t_eval, inputs={"rate": 0.1}, calculate_sensitivities=True + ) + assert len(solution_sens.sensitivities) == 0