#247 merge master

pybamm-team · Feb 21, 2020 · ccae1df · ccae1df
2 parents 36aa13b + c1101ba
commit ccae1df
Show file tree

Hide file tree

Showing 11 changed files with 238 additions and 120 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -2,6 +2,7 @@
 
 ## Features
 
+-   Changed rootfinding algorithm to CasADi, scipy.optimize.root still accessible as an option ([#844](https://github.com/pybamm-team/PyBaMM/pull/844))
 -   Added capacitance effects to lithium-ion models ([#842](https://github.com/pybamm-team/PyBaMM/pull/842))
 -   Added NCA parameter set ([#824](https://github.com/pybamm-team/PyBaMM/pull/824))
 -   Added functionality to `Solution` that automatically gets `t_eval` from the data when simulating drive cycles and performs checks to ensure the output has the required resolution to accurately capture the input current ([#819](https://github.com/pybamm-team/PyBaMM/pull/819))

diff --git a/pybamm/models/full_battery_models/lithium_ion/basic_dfn.py b/pybamm/models/full_battery_models/lithium_ion/basic_dfn.py
@@ -177,11 +177,17 @@ def __init__(self, name="Doyle-Fuller-Newman model"):
         # Boundary conditions must be provided for equations with spatial derivatives
         self.boundary_conditions[c_s_n] = {
             "left": (pybamm.Scalar(0), "Neumann"),
-            "right": (-param.C_n * j_n / param.a_n, "Neumann"),
+            "right": (
+                -param.C_n * j_n / param.a_n / param.D_n(c_s_surf_n, T),
+                "Neumann",
+            ),
         }
         self.boundary_conditions[c_s_p] = {
             "left": (pybamm.Scalar(0), "Neumann"),
-            "right": (-param.C_p * j_p / param.a_p / param.gamma_p, "Neumann"),
+            "right": (
+                -param.C_p * j_p / param.a_p / param.gamma_p / param.D_p(c_s_surf_p, T),
+                "Neumann",
+            ),
         }
         # c_n_init and c_p_init can in general be functions of x
         # Note the broadcasting, for domains
@@ -195,14 +201,22 @@ def __init__(self, name="Doyle-Fuller-Newman model"):
         self.initial_conditions[c_s_p] = param.c_p_init(x_p)
         # Events specify points at which a solution should terminate
         self.events += [
-            pybamm.Event("Minimum negative particle surface concentration",
-                         pybamm.min(c_s_surf_n) - 0.01),
-            pybamm.Event("Maximum negative particle surface concentration",
-                         (1 - 0.01) - pybamm.max(c_s_surf_n)),
-            pybamm.Event("Minimum positive particle surface concentration",
-                         pybamm.min(c_s_surf_p) - 0.01),
-            pybamm.Event("Maximum positive particle surface concentration",
-                         (1 - 0.01) - pybamm.max(c_s_surf_p)),
+            pybamm.Event(
+                "Minimum negative particle surface concentration",
+                pybamm.min(c_s_surf_n) - 0.01,
+            ),
+            pybamm.Event(
+                "Maximum negative particle surface concentration",
+                (1 - 0.01) - pybamm.max(c_s_surf_n),
+            ),
+            pybamm.Event(
+                "Minimum positive particle surface concentration",
+                pybamm.min(c_s_surf_p) - 0.01,
+            ),
+            pybamm.Event(
+                "Maximum positive particle surface concentration",
+                (1 - 0.01) - pybamm.max(c_s_surf_p),
+            ),
         ]
         ######################
         # Current in the solid
@@ -257,8 +271,11 @@ def __init__(self, name="Doyle-Fuller-Newman model"):
             "right": (pybamm.Scalar(0), "Neumann"),
         }
         self.initial_conditions[c_e] = param.c_e_init
-        self.events.append(pybamm.Event("Zero electrolyte concentration cut-off",
-                                        pybamm.min(c_e) - 0.002))
+        self.events.append(
+            pybamm.Event(
+                "Zero electrolyte concentration cut-off", pybamm.min(c_e) - 0.002
+            )
+        )
 
         ######################
         # (Some) variables

diff --git a/pybamm/models/full_battery_models/lithium_ion/basic_spm.py b/pybamm/models/full_battery_models/lithium_ion/basic_spm.py
@@ -80,34 +80,48 @@ def __init__(self, name="Single Particle Model"):
         N_s_p = -param.D_p(c_s_p, T) * pybamm.grad(c_s_p)
         self.rhs[c_s_n] = -(1 / param.C_n) * pybamm.div(N_s_n)
         self.rhs[c_s_p] = -(1 / param.C_p) * pybamm.div(N_s_p)
+        # Surf takes the surface value of a variable, i.e. its boundary value on the
+        # right side. This is also accessible via `boundary_value(x, "right")`, with
+        # "left" providing the boundary value of the left side
+        c_s_surf_n = pybamm.surf(c_s_n)
+        c_s_surf_p = pybamm.surf(c_s_p)
         # Boundary conditions must be provided for equations with spatial derivatives
         self.boundary_conditions[c_s_n] = {
             "left": (pybamm.Scalar(0), "Neumann"),
-            "right": (-param.C_n * j_n / param.a_n, "Neumann"),
+            "right": (
+                -param.C_n * j_n / param.a_n / param.D_n(c_s_surf_n, T),
+                "Neumann",
+            ),
         }
         self.boundary_conditions[c_s_p] = {
             "left": (pybamm.Scalar(0), "Neumann"),
-            "right": (-param.C_p * j_p / param.a_p / param.gamma_p, "Neumann"),
+            "right": (
+                -param.C_p * j_p / param.a_p / param.gamma_p / param.D_p(c_s_surf_p, T),
+                "Neumann",
+            ),
         }
         # c_n_init and c_p_init are functions, but for the SPM we evaluate them at x=0
         # and x=1 since there is no x-dependence in the particles
         self.initial_conditions[c_s_n] = param.c_n_init(0)
         self.initial_conditions[c_s_p] = param.c_p_init(1)
-        # Surf takes the surface value of a variable, i.e. its boundary value on the
-        # right side. This is also accessible via `boundary_value(x, "right")`, with
-        # "left" providing the boundary value of the left side
-        c_s_surf_n = pybamm.surf(c_s_n)
-        c_s_surf_p = pybamm.surf(c_s_p)
         # Events specify points at which a solution should terminate
         self.events += [
-            pybamm.Event("Minimum negative particle surface concentration",
-                         pybamm.min(c_s_surf_n) - 0.01),
-            pybamm.Event("Maximum negative particle surface concentration",
-                         (1 - 0.01) - pybamm.max(c_s_surf_n)),
-            pybamm.Event("Minimum positive particle surface concentration",
-                         pybamm.min(c_s_surf_p) - 0.01),
-            pybamm.Event("Maximum positive particle surface concentration",
-                         (1 - 0.01) - pybamm.max(c_s_surf_p)),
+            pybamm.Event(
+                "Minimum negative particle surface concentration",
+                pybamm.min(c_s_surf_n) - 0.01,
+            ),
+            pybamm.Event(
+                "Maximum negative particle surface concentration",
+                (1 - 0.01) - pybamm.max(c_s_surf_n),
+            ),
+            pybamm.Event(
+                "Minimum positive particle surface concentration",
+                pybamm.min(c_s_surf_p) - 0.01,
+            ),
+            pybamm.Event(
+                "Maximum positive particle surface concentration",
+                (1 - 0.01) - pybamm.max(c_s_surf_p),
+            ),
         ]
 
         # Note that the SPM does not have any algebraic equations, so the `algebraic`