Merge branch 'develop' into remove-unneeded-params

Author: valentinsulzer

CHANGELOG.md

@@ -9,6 +9,7 @@
 
 ## Bug fixes
 
+- Fixed electrolyte conservation in the case of concentration-dependent transference number ([#2758](https://github.com/pybamm-team/PyBaMM/pull/2758))
 - Fixed `plot_voltage_components` so that the sum of overpotentials is now equal to the voltage ([#2740](https://github.com/pybamm-team/PyBaMM/pull/2740))
 
 ## Optimizations
@@ -18,6 +19,7 @@
 ## Breaking changes
 
 - Removed "... cation signed stoichiometry" and "... electrons in reaction" parameters, they are now hardcoded. ([#2778](https://github.com/pybamm-team/PyBaMM/pull/2778))
+- When using `solver.step()`, the first time point in the step is shifted by `pybamm.settings.step_start_offset` (default 1 ns) to avoid having duplicate times in the solution steps from the end of one step and the start of the next. ([#2773](https://github.com/pybamm-team/PyBaMM/pull/2773))
 - Renamed "Measured open circuit voltage [V]" to "Surface open-circuit voltage [V]". This variable was calculated from surface particle concentrations, and hence "hid" the overpotential from particle gradients. The new variable "Open-circuit voltage [V]" is calculated from bulk particle concentrations instead. ([#2740](https://github.com/pybamm-team/PyBaMM/pull/2740))
 - Renamed all references to "open circuit" to be "open-circuit" instead. ([#2740](https://github.com/pybamm-team/PyBaMM/pull/2740))
 - Renamed parameter "1 + dlnf/dlnc" to "Thermodynamic factor". ([#2727](https://github.com/pybamm-team/PyBaMM/pull/2727))
@@ -30,6 +32,7 @@
 - Added an option for using a banded jacobian and sundials banded solvers for the IDAKLU solve ([#2677](https://github.com/pybamm-team/PyBaMM/pull/2677))
 - The "particle size" option can now be a tuple to allow different behaviour in each electrode ([#2672](https://github.com/pybamm-team/PyBaMM/pull/2672)).
 - Added temperature control to experiment class. ([#2518](https://github.com/pybamm-team/PyBaMM/pull/2518))
+- Added method to calculate maximum theoretical energy. ([#2777](https://github.com/pybamm-team/PyBaMM/pull/2777))
 
 ## Bug fixes
 

examples/notebooks/Getting Started/Tutorial 7 - Model options.ipynb

@@ -59,7 +59,7 @@
     {
      "data": {
       "text/plain": [
-       "<pybamm.solvers.solution.Solution at 0x7f29f1333e50>"
+       "<pybamm.solvers.solution.Solution at 0x139461520>"
       ]
      },
      "execution_count": 3,
@@ -70,10 +70,7 @@
    "source": [
     "model = pybamm.lithium_ion.SPMe(options=options) # loading in options\n",
     "\n",
-    "parameter_values = model.default_parameter_values\n",
-    "parameter_values[\"Current function [A]\"] = 3\n",
-    "\n",
-    "sim = pybamm.Simulation(model, parameter_values=parameter_values)\n",
+    "sim = pybamm.Simulation(model)\n",
     "sim.solve([0, 3600])"
    ]
   },
@@ -92,16 +89,26 @@
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "c00d126a3fc543829596dd24f8b9d626",
+       "model_id": "0fc5db759e804d3fb2793222ddbec5d4",
        "version_major": 2,
        "version_minor": 0
       },
       "text/plain": [
-       "interactive(children=(FloatSlider(value=0.0, description='t', max=800.0000000000001, step=8.000000000000002), …"
+       "interactive(children=(FloatSlider(value=0.0, description='t', max=1.0, step=0.01), Output()), _dom_classes=('w…"
       ]
      },
      "metadata": {},
      "output_type": "display_data"
+    },
+    {
+     "data": {
+      "text/plain": [
+       "<pybamm.plotting.quick_plot.QuickPlot at 0x139461fd0>"
+      ]
+     },
+     "execution_count": 4,
+     "metadata": {},
+     "output_type": "execute_result"
     }
    ],
    "source": [
@@ -139,7 +146,8 @@
       "[1] Joel A. E. Andersson, Joris Gillis, Greg Horn, James B. Rawlings, and Moritz Diehl. CasADi – A software framework for nonlinear optimization and optimal control. Mathematical Programming Computation, 11(1):1–36, 2019. doi:10.1007/s12532-018-0139-4.\n",
       "[2] Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, and others. Array programming with NumPy. Nature, 585(7825):357–362, 2020. doi:10.1038/s41586-020-2649-2.\n",
       "[3] Scott G. Marquis, Valentin Sulzer, Robert Timms, Colin P. Please, and S. Jon Chapman. An asymptotic derivation of a single particle model with electrolyte. Journal of The Electrochemical Society, 166(15):A3693–A3706, 2019. doi:10.1149/2.0341915jes.\n",
-      "[4] Valentin Sulzer, Scott G. Marquis, Robert Timms, Martin Robinson, and S. Jon Chapman. Python Battery Mathematical Modelling (PyBaMM). ECSarXiv. February, 2020. doi:10.1149/osf.io/67ckj.\n",
+      "[4] Valentin Sulzer, Scott G. Marquis, Robert Timms, Martin Robinson, and S. Jon Chapman. Python Battery Mathematical Modelling (PyBaMM). Journal of Open Research Software, 9(1):14, 2021. doi:10.5334/jors.309.\n",
+      "[5] Robert Timms, Scott G Marquis, Valentin Sulzer, Colin P. Please, and S Jonathan Chapman. Asymptotic Reduction of a Lithium-ion Pouch Cell Model. SIAM Journal on Applied Mathematics, 81(3):765–788, 2021. doi:10.1137/20M1336898.\n",
       "\n"
      ]
     }
@@ -165,7 +173,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.1"
+   "version": "3.9.16"
   },
   "vscode": {
    "interpreter": {

pybamm/models/full_battery_models/lithium_ion/electrode_soh.py

@@ -557,3 +557,47 @@ def get_min_max_stoichiometries(
     """
     esoh_solver = ElectrodeSOHSolver(parameter_values, param, known_value)
     return esoh_solver.get_min_max_stoichiometries()
+
+
+def calculate_theoretical_energy(
+    parameter_values, initial_soc=1.0, final_soc=0.0, points=100
+):
+    """
+    Calculate maximum energy possible from a cell given OCV, initial soc, and final soc
+    given voltage limits, open-circuit potentials, etc defined by parameter_values
+
+    Parameters
+    ----------
+    parameter_values : :class:`pybamm.ParameterValues`
+        The parameter values class that will be used for the simulation.
+    initial_soc : float
+        The soc at begining of discharge, default 1.0
+    final_soc : float
+        The soc at end of discharge, default 1.0
+    points : int
+        The number of points at which to calculate voltage.
+
+    Returns
+    -------
+    E
+        The total energy of the cell in Wh
+    """
+    # Get initial and final stoichiometric values.
+    n_i, p_i = get_initial_stoichiometries(initial_soc, parameter_values)
+    n_f, p_f = get_initial_stoichiometries(final_soc, parameter_values)
+    n_vals = np.linspace(n_i, n_f, num=points)
+    p_vals = np.linspace(p_i, p_f, num=points)
+    # Calculate OCV at each stoichiometry
+    param = pybamm.LithiumIonParameters()
+    T = param.T_amb(0)
+    Vs = np.empty(n_vals.shape)
+    for i in range(n_vals.size):
+        Vs[i] = parameter_values.evaluate(
+            param.p.prim.U(p_vals[i], T)
+        ) - parameter_values.evaluate(param.n.prim.U(n_vals[i], T))
+    # Calculate dQ
+    Q_p = parameter_values.evaluate(param.p.prim.Q_init) * (p_f - p_i)
+    dQ = Q_p / (points - 1)
+    # Integrate and convert to W-h
+    E = np.trapz(Vs, dx=dQ)
+    return E

pybamm/models/submodels/electrolyte_diffusion/full_diffusion.py

@@ -89,17 +89,12 @@ def get_coupled_variables(self, variables):
     def set_rhs(self, variables):
         eps_c_e = variables["Porosity times concentration [mol.m-3]"]
         c_e = variables["Electrolyte concentration [mol.m-3]"]
-        T = variables["Cell temperature [K]"]
-        N_e_diffusion = variables["Electrolyte diffusion flux [mol.m-2.s-1]"]
-        N_e_convection = variables["Electrolyte convection flux [mol.m-2.s-1]"]
+        N_e = variables["Electrolyte flux [mol.m-2.s-1]"]
         div_Vbox = variables["Transverse volume-averaged acceleration [m.s-2]"]
 
         sum_s_a_j = variables["Sum of electrolyte reaction source terms [A.m-3]"]
-        sum_a_j = variables["Sum of volumetric interfacial current densities [A.m-3]"]
         sum_s_a_j.print_name = "aj"
-        source_terms = (sum_s_a_j - self.param.t_plus(c_e, T) * sum_a_j) / self.param.F
-
-        N_e = N_e_diffusion + N_e_convection
+        source_terms = sum_s_a_j / self.param.F
 
         self.rhs = {eps_c_e: -pybamm.div(N_e) + source_terms - c_e * div_Vbox}
 

pybamm/settings.py

@@ -12,6 +12,7 @@ class Settings(object):
     _abs_smoothing = "exact"
     max_words_in_line = 4
     max_y_value = 1e5
+    step_start_offset = 1e-9
     tolerances = {
         "D_e__c_e": 10,  # dimensional
         "kappa_e__c_e": 10,  # dimensional

pybamm/solvers/base_solver.py

@@ -1085,18 +1085,34 @@ def step(
                     "Cannot step empty model, use `pybamm.DummySolver` instead"
                 )
 
-        # Make sure dt is positive
-        if dt <= 0:
-            raise pybamm.SolverError("Step time must be positive")
+        # Make sure dt is greater than the offset
+        step_start_offset = pybamm.settings.step_start_offset
+        if dt <= step_start_offset:
+            raise pybamm.SolverError(
+                f"Step time must be at least {pybamm.TimerTime(step_start_offset)}"
+            )
+
+        t_start = old_solution.t[-1]
+        t_end = t_start + dt
+        # Calculate t_eval
+        t_eval = np.linspace(t_start, t_end, npts)
+
+        if t_start == 0:
+            t_start_shifted = t_start
+        else:
+            # offset t_start by t_start_offset (default 1 ns)
+            # to avoid repeated times in the solution
+            # from having the same time at the end of the previous step and
+            # the start of the next step
+            t_start_shifted = t_start + step_start_offset
+            t_eval[0] = t_start_shifted
 
         # Set timer
         timer = pybamm.Timer()
 
         # Set up inputs
         model_inputs = self._set_up_model_inputs(model, inputs)
 
-        t = old_solution.t[-1]
-
         first_step_this_model = False
         if model not in self._model_set_up:
             first_step_this_model = True
@@ -1139,16 +1155,17 @@ def step(
         set_up_time = timer.time()
 
         # (Re-)calculate consistent initial conditions
-        self._set_initial_conditions(model, t, model_inputs, update_rhs=False)
-
-        # Calculate t_eval
-        t_eval = np.linspace(t, t + dt, npts)
+        self._set_initial_conditions(
+            model, t_start_shifted, model_inputs, update_rhs=False
+        )
 
         # Check initial conditions don't violate events
         self._check_events_with_initial_conditions(t_eval, model, model_inputs)
 
         # Step
-        pybamm.logger.verbose("Stepping for {:.0f} < t < {:.0f}".format(t, (t + dt)))
+        pybamm.logger.verbose(
+            "Stepping for {:.0f} < t < {:.0f}".format(t_start_shifted, t_end)
+        )
         timer.reset()
         solution = self._integrate(model, t_eval, model_inputs)
         solution.solve_time = timer.time()

tests/integration/test_experiments.py

@@ -14,20 +14,21 @@ def test_discharge_rest_charge(self):
                 "Rest for 1 hour",
                 "Charge at C/2 for 1 hour",
             ],
-            period="0.25 hours",
+            period="0.5 hours",
         )
         model = pybamm.lithium_ion.SPM()
         sim = pybamm.Simulation(
             model, experiment=experiment, solver=pybamm.CasadiSolver()
         )
         sim.solve()
         np.testing.assert_array_almost_equal(
-            sim._solution["Time [h]"].entries, np.linspace(0, 3, 13)
+            sim._solution["Time [h]"].entries,
+            np.array([0, 0.5, 1, 1 + 1e-9, 1.5, 2, 2 + 1e-9, 2.5, 3]),
         )
         cap = model.default_parameter_values["Nominal cell capacity [A.h]"]
         np.testing.assert_array_almost_equal(
             sim._solution["Current [A]"].entries,
-            [cap / 2] * 5 + [0] * 4 + [-cap / 2] * 4,
+            [cap / 2] * 3 + [0] * 3 + [-cap / 2] * 3,
         )
 
     def test_rest_discharge_rest(self):
@@ -57,13 +58,10 @@ def test_gitt(self):
             model, experiment=experiment, solver=pybamm.CasadiSolver()
         )
         sim.solve()
-        np.testing.assert_array_almost_equal(
-            sim._solution["Time [h]"].entries, np.arange(0, 20.01, 0.1)
-        )
         cap = model.default_parameter_values["Nominal cell capacity [A.h]"]
         np.testing.assert_array_almost_equal(
             sim._solution["Current [A]"].entries,
-            [cap / 20] * 11 + [0] * 10 + ([cap / 20] * 10 + [0] * 10) * 9,
+            [cap / 20] * 11 + [0] * 11 + ([cap / 20] * 11 + [0] * 11) * 9,
         )
 
     def test_infeasible(self):

tests/integration/test_models/standard_output_comparison.py

@@ -128,13 +128,13 @@ def test_all(self):
         self.compare("X-averaged positive electrode open-circuit potential [V]")
         self.compare("Voltage [V]")
         self.compare("X-averaged solid phase ohmic losses [V]")
-        self.compare("Negative electrode reaction overpotential [V]")
+        self.compare("Negative electrode reaction overpotential [V]", atol=1e-4)
         self.compare("Positive electrode reaction overpotential [V]")
         self.compare("Negative electrode potential [V]", atol=1e-5)
         self.compare("Positive electrode potential [V]")
         self.compare("Electrolyte potential [V]")
         # Currents
-        self.compare("Exchange current density [A.m-2]")
+        self.compare("Exchange current density [A.m-2]", atol=2e-3)
         self.compare("Negative electrode current density [A.m-2]", atol=1e-10)
         self.compare("Positive electrode current density [A.m-2]", atol=1e-10)
 

tests/integration/test_models/standard_output_tests.py

@@ -545,7 +545,10 @@ def test_conservation(self):
         diff = (
             self.c_e_tot(self.solution.t[1:]) - self.c_e_tot(self.solution.t[:-1])
         ) / self.c_e_tot(self.solution.t[:-1])
-        if self.model.options["surface form"] == "differential":
+        if self.model.options["surface form"] == "differential" or (
+            isinstance(self.model, pybamm.lithium_ion.DFN)
+            and self.model.options["surface form"] == "algebraic"
+        ):
             np.testing.assert_allclose(0, diff, atol=1e-4, rtol=1e-6)
         else:
             np.testing.assert_allclose(0, diff, atol=1e-14, rtol=1e-14)

tests/unit/test_experiments/test_simulation_with_experiment.py

@@ -176,22 +176,18 @@ def test_run_experiment_cccv_ode(self):
             solution = sim.solve(solver=pybamm.CasadiSolver("fast with events"))
             solutions.append(solution)
 
+        t = solutions[1]["Time [s]"].data
         np.testing.assert_array_almost_equal(
-            solutions[0]["Voltage [V]"].data,
-            solutions[1]["Voltage [V]"].data,
+            solutions[0]["Voltage [V]"](t=t),
+            solutions[1]["Voltage [V]"](t=t),
             decimal=1,
         )
         np.testing.assert_array_almost_equal(
-            solutions[0]["Current [A]"].data,
-            solutions[1]["Current [A]"].data,
+            solutions[0]["Current [A]"](t=t),
+            solutions[1]["Current [A]"](t=t),
             decimal=0,
         )
 
-        np.testing.assert_array_equal(
-            solutions[0]["Ambient temperature [C]"].data,
-            solutions[1]["Ambient temperature [C]"].data,
-        )
-
         self.assertEqual(solutions[1].termination, "final time")
 
     @unittest.skipIf(not pybamm.have_idaklu(), "idaklu solver is not installed")
@@ -516,12 +512,13 @@ def test_run_experiment_skip_steps(self):
         )
         sol2 = sim2.solve()
         np.testing.assert_array_almost_equal(
-            sol["Voltage [V]"].data, sol2["Voltage [V]"].data
+            sol["Voltage [V]"].data, sol2["Voltage [V]"].data, decimal=5
         )
         for idx1, idx2 in [(1, 0), (2, 1), (4, 2)]:
             np.testing.assert_array_almost_equal(
                 sol.cycles[0].steps[idx1]["Voltage [V]"].data,
                 sol2.cycles[0].steps[idx2]["Voltage [V]"].data,
+                decimal=5,
             )
 
     def test_all_empty_solution_errors(self):