Merge pull request #786 from pybamm-team/issue-785-init-concs

Issue 785 init concs

CHANGELOG.md

@@ -5,6 +5,7 @@
 -   Added fuzzy string matching for parameters and variables ([#796](https://github.com/pybamm-team/PyBaMM/pull/796))
 -   Changed ParameterValues to raise an error when a parameter that wasn't previously defined is updated ([#796](https://github.com/pybamm-team/PyBaMM/pull/796))
 -   Added some basic models (BasicSPM and BasicDFN) in order to clearly demonstrate the PyBaMM model structure for battery models ([#795](https://github.com/pybamm-team/PyBaMM/pull/795))
+-   Allow initial conditions in the particle to depend on x ([#786](https://github.com/pybamm-team/PyBaMM/pull/786))
 -   Added the harmonic mean to the Finite Volume method, which is now used when computing fluxes ([#783](https://github.com/pybamm-team/PyBaMM/pull/783))
 -   Refactored `Solution` to make it a dictionary that contains all of the solution variables. This automatically creates `ProcessedVariable` objects when required, so that the solution can be obtained much more easily. ([#781](https://github.com/pybamm-team/PyBaMM/pull/781))
 -   Added notebook to explain broadcasts ([#776](https://github.com/pybamm-team/PyBaMM/pull/776))

pybamm/expression_tree/parameter.py

@@ -1,6 +1,7 @@
 #
 # Parameter classes
 #
+import numbers
 import numpy as np
 import pybamm
 
@@ -60,10 +61,19 @@ def __init__(self, name, *children, diff_variable=None):
         # assign diff variable
         self.diff_variable = diff_variable
         children_list = list(children)
+
+        # Turn numbers into scalars
+        for idx, child in enumerate(children_list):
+            if isinstance(child, numbers.Number):
+                children_list[idx] = pybamm.Scalar(child)
+
         domain = self.get_children_domains(children_list)
-        auxiliary_domains = self.get_children_auxiliary_domains(children)
+        auxiliary_domains = self.get_children_auxiliary_domains(children_list)
         super().__init__(
-            name, children=children, domain=domain, auxiliary_domains=auxiliary_domains
+            name,
+            children=children_list,
+            domain=domain,
+            auxiliary_domains=auxiliary_domains,
         )
 
     def set_id(self):

pybamm/models/full_battery_models/lithium_ion/basic_dfn.py

@@ -177,8 +177,16 @@ def __init__(self, name="Doyle-Fuller-Newman model"):
             "left": (pybamm.Scalar(0), "Neumann"),
             "right": (-param.C_p * j_p / param.a_p / param.gamma_p, "Neumann"),
         }
-        self.initial_conditions[c_s_n] = param.c_n_init
-        self.initial_conditions[c_s_p] = param.c_p_init
+        # c_n_init and c_p_init can in general be functions of x
+        # Note the broadcasting, for domains
+        x_n = pybamm.PrimaryBroadcast(
+            pybamm.standard_spatial_vars.x_n, "negative particle"
+        )
+        self.initial_conditions[c_s_n] = param.c_n_init(x_n)
+        x_p = pybamm.PrimaryBroadcast(
+            pybamm.standard_spatial_vars.x_p, "positive particle"
+        )
+        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",
@@ -211,10 +219,12 @@ def __init__(self, name="Doyle-Fuller-Newman model"):
         # Initial conditions must also be provided for algebraic equations, as an
         # initial guess for a root-finding algorithm which calculates consistent initial
         # conditions
+        # We evaluate c_n_init at x=0 and c_p_init at x=1 (this is just an initial
+        # guess so actual value is not too important)
         self.initial_conditions[phi_s_n] = pybamm.Scalar(0)
         self.initial_conditions[phi_s_p] = param.U_p(
-            param.c_p_init, param.T_init
-        ) - param.U_n(param.c_n_init, param.T_init)
+            param.c_p_init(1), param.T_init
+        ) - param.U_n(param.c_n_init(0), param.T_init)
 
         ######################
         # Current in the electrolyte
@@ -227,7 +237,7 @@ def __init__(self, name="Doyle-Fuller-Newman model"):
             "left": (pybamm.Scalar(0), "Neumann"),
             "right": (pybamm.Scalar(0), "Neumann"),
         }
-        self.initial_conditions[phi_e] = -param.U_n(param.c_n_init, param.T_init)
+        self.initial_conditions[phi_e] = -param.U_n(param.c_n_init(0), param.T_init)
 
         ######################
         # Electrolyte concentration

pybamm/models/full_battery_models/lithium_ion/basic_spm.py

@@ -89,8 +89,10 @@ def __init__(self, name="Single Particle Model"):
             "left": (pybamm.Scalar(0), "Neumann"),
             "right": (-param.C_p * j_p / param.a_p / param.gamma_p, "Neumann"),
         }
-        self.initial_conditions[c_s_n] = param.c_n_init
-        self.initial_conditions[c_s_p] = param.c_p_init
+        # 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

pybamm/models/submodels/electrode/ohm/full_ohm.py

@@ -95,8 +95,8 @@ def set_initial_conditions(self, variables):
         if self.domain == "Negative":
             phi_s_init = pybamm.Scalar(0)
         elif self.domain == "Positive":
-            phi_s_init = self.param.U_p(self.param.c_p_init, T_init) - self.param.U_n(
-                self.param.c_n_init, T_init
-            )
+            phi_s_init = self.param.U_p(
+                self.param.c_p_init(1), T_init
+            ) - self.param.U_n(self.param.c_n_init(0), T_init)
 
         self.initial_conditions[phi_s] = phi_s_init

pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/full_stefan_maxwell_conductivity.py

@@ -64,4 +64,6 @@ def set_algebraic(self, variables):
     def set_initial_conditions(self, variables):
         phi_e = variables["Electrolyte potential"]
         T_init = self.param.T_init
-        self.initial_conditions = {phi_e: -self.param.U_n(self.param.c_n_init, T_init)}
+        self.initial_conditions = {
+            phi_e: -self.param.U_n(self.param.c_n_init(0), T_init)
+        }

pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/surface_potential_form/full_surface_form_stefan_maxwell_conductivity.py

@@ -45,9 +45,9 @@ def set_initial_conditions(self, variables):
 
         delta_phi_e = variables[self.domain + " electrode surface potential difference"]
         if self.domain == "Negative":
-            delta_phi_e_init = self.param.U_n(self.param.c_n_init, self.param.T_init)
+            delta_phi_e_init = self.param.U_n(self.param.c_n_init(0), self.param.T_init)
         elif self.domain == "Positive":
-            delta_phi_e_init = self.param.U_p(self.param.c_p_init, self.param.T_init)
+            delta_phi_e_init = self.param.U_p(self.param.c_p_init(1), self.param.T_init)
 
         self.initial_conditions = {delta_phi_e: delta_phi_e_init}
 

pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/surface_potential_form/leading_surface_form_stefan_maxwell_conductivity.py

@@ -49,9 +49,9 @@ def set_initial_conditions(self, variables):
             + " electrode surface potential difference"
         ]
         if self.domain == "Negative":
-            delta_phi_init = self.param.U_n(self.param.c_n_init, self.param.T_init)
+            delta_phi_init = self.param.U_n(self.param.c_n_init(0), self.param.T_init)
         elif self.domain == "Positive":
-            delta_phi_init = self.param.U_p(self.param.c_p_init, self.param.T_init)
+            delta_phi_init = self.param.U_p(self.param.c_p_init(1), self.param.T_init)
 
         self.initial_conditions = {delta_phi: delta_phi_init}
 

pybamm/models/submodels/particle/base_particle.py

@@ -76,17 +76,6 @@ def _flux_law(self, c, T):
     def _unpack(self, variables):
         raise NotImplementedError
 
-    def set_initial_conditions(self, variables):
-        c, _, _ = self._unpack(variables)
-
-        if self.domain == "Negative":
-            c_init = self.param.c_n_init
-
-        elif self.domain == "Positive":
-            c_init = self.param.c_p_init
-
-        self.initial_conditions = {c: c_init}
-
     def set_events(self, variables):
         c_s_surf = variables[self.domain + " particle surface concentration"]
         tol = 0.01

pybamm/models/submodels/particle/fast/base_fast_particle.py

@@ -2,7 +2,6 @@
 # Base class for particles each with uniform concentration (i.e. infinitely fast
 # diffusion in r)
 #
-
 from ..base_particle import BaseParticle
 
 

pybamm/models/submodels/particle/fast/fast_many_particles.py

@@ -69,3 +69,16 @@ def _unpack(self, variables):
         j = variables[self.domain + " electrode interfacial current density"]
 
         return c_s_surf, N_s, j
+
+    def set_initial_conditions(self, variables):
+        c, _, _ = self._unpack(variables)
+
+        if self.domain == "Negative":
+            x_n = pybamm.standard_spatial_vars.x_n
+            c_init = self.param.c_n_init(x_n)
+
+        elif self.domain == "Positive":
+            x_p = pybamm.standard_spatial_vars.x_p
+            c_init = self.param.c_p_init(x_p)
+
+        self.initial_conditions = {c: c_init}

pybamm/models/submodels/particle/fast/fast_single_particle.py

@@ -77,3 +77,19 @@ def _unpack(self, variables):
         ]
 
         return c_s_surf_xav, N_s_xav, j_av
+
+    def set_initial_conditions(self, variables):
+        """
+        For single particle models, initial conditions can't depend on x so we
+        arbitrarily evaluate them at x=0 in the negative electrode and x=1 in the
+        positive electrode (they will usually be constant)
+        """
+        c, _, _ = self._unpack(variables)
+
+        if self.domain == "Negative":
+            c_init = self.param.c_n_init(0)
+
+        elif self.domain == "Positive":
+            c_init = self.param.c_p_init(1)
+
+        self.initial_conditions = {c: c_init}

pybamm/models/submodels/particle/fickian/fickian_many_particles.py

@@ -76,3 +76,20 @@ def _unpack(self, variables):
         j = variables[self.domain + " electrode interfacial current density"]
 
         return c_s, N_s, j
+
+    def set_initial_conditions(self, variables):
+        c, _, _ = self._unpack(variables)
+
+        if self.domain == "Negative":
+            x_n = pybamm.PrimaryBroadcast(
+                pybamm.standard_spatial_vars.x_n, "negative particle"
+            )
+            c_init = self.param.c_n_init(x_n)
+
+        elif self.domain == "Positive":
+            x_p = pybamm.PrimaryBroadcast(
+                pybamm.standard_spatial_vars.x_p, "positive particle"
+            )
+            c_init = self.param.c_p_init(x_p)
+
+        self.initial_conditions = {c: c_init}

pybamm/models/submodels/particle/fickian/fickian_single_particle.py

@@ -74,3 +74,19 @@ def _unpack(self, variables):
         ]
 
         return c_s_xav, N_s_xav, j_av
+
+    def set_initial_conditions(self, variables):
+        """
+        For single particle models, initial conditions can't depend on x so we
+        arbitrarily evaluate them at x=0 in the negative electrode and x=1 in the
+        positive electrode (they will usually be constant)
+        """
+        c, _, _ = self._unpack(variables)
+
+        if self.domain == "Negative":
+            c_init = self.param.c_n_init(0)
+
+        elif self.domain == "Positive":
+            c_init = self.param.c_p_init(1)
+
+        self.initial_conditions = {c: c_init}

pybamm/parameters/standard_parameters_lead_acid.py

@@ -420,9 +420,13 @@ def U_p_dimensional(c_e, T):
 
 
 # hack to make consistent ic with lithium-ion
-# find a way to not have to do this
-c_n_init = c_e_init
-c_p_init = c_e_init
+def c_n_init(x):
+    return c_e_init
+
+
+def c_p_init(x):
+    return c_e_init
+
 
 # Thermal effects not implemented for lead-acid, but parameters needed for consistency
 T_init = pybamm.Scalar(0)

pybamm/parameters/standard_parameters_lithium_ion.py

@@ -99,12 +99,21 @@
 c_e_init_dimensional = pybamm.Parameter(
     "Initial concentration in electrolyte [mol.m-3]"
 )
-c_n_init_dimensional = pybamm.Parameter(
-    "Initial concentration in negative electrode [mol.m-3]"
-)
-c_p_init_dimensional = pybamm.Parameter(
-    "Initial concentration in positive electrode [mol.m-3]"
-)
+
+
+def c_n_init_dimensional(x):
+    "Initial concentration as a function of dimensionless position x"
+    return pybamm.FunctionParameter(
+        "Initial concentration in negative electrode [mol.m-3]", x
+    )
+
+
+def c_p_init_dimensional(x):
+    "Initial concentration as a function of dimensionless position x"
+    return pybamm.FunctionParameter(
+        "Initial concentration in positive electrode [mol.m-3]", x
+    )
+
 
 # thermal
 Delta_T = pybamm.thermal_parameters.Delta_T
@@ -362,10 +371,18 @@ def chi(c_e):
 )
 
 # Initial conditions
-c_e_init = c_e_init_dimensional / c_e_typ
-c_n_init = c_n_init_dimensional / c_n_max
-c_p_init = c_p_init_dimensional / c_p_max
 T_init = pybamm.thermal_parameters.T_init
+c_e_init = c_e_init_dimensional / c_e_typ
+
+
+def c_n_init(x):
+    "Dimensionless initial concentration as a function of dimensionless position x"
+    return c_n_init_dimensional(x) / c_n_max
+
+
+def c_p_init(x):
+    "Dimensionless initial concentration as a function of dimensionless position x"
+    return c_p_init_dimensional(x) / c_p_max
 
 
 # --------------------------------------------------------------------------------------

tests/unit/test_models/test_submodels/test_particle/test_base_particle.py

@@ -9,15 +9,17 @@
 
 class TestBaseParticle(unittest.TestCase):
     def test_public_functions(self):
+        variables = {
+            "Negative particle surface concentration": 0,
+            "Positive particle surface concentration": 0,
+        }
         submodel = pybamm.particle.BaseParticle(None, "Negative")
-        std_tests = tests.StandardSubModelTests(submodel)
-        with self.assertRaises(NotImplementedError):
-            std_tests.test_all()
+        std_tests = tests.StandardSubModelTests(submodel, variables)
+        std_tests.test_all()
 
         submodel = pybamm.particle.BaseParticle(None, "Positive")
-        std_tests = tests.StandardSubModelTests(submodel)
-        with self.assertRaises(NotImplementedError):
-            std_tests.test_all()
+        std_tests = tests.StandardSubModelTests(submodel, variables)
+        std_tests.test_all()
 
 
 if __name__ == "__main__":

tests/unit/test_parameters/test_dimensionless_parameter_values_lithium_ion.py

@@ -63,7 +63,7 @@ def test_lithium_ion(self):
         "particle dynamics"
         # neg diffusion coefficient
         np.testing.assert_almost_equal(
-            values.evaluate(param.D_n_dimensional(param.c_n_init, param.T_ref)),
+            values.evaluate(param.D_n_dimensional(param.c_n_init(0), param.T_ref)),
             3.9 * 10 ** (-14),
             2,
         )
@@ -78,7 +78,7 @@ def test_lithium_ion(self):
 
         # pos diffusion coefficient
         np.testing.assert_almost_equal(
-            values.evaluate(param.D_p_dimensional(param.c_p_init, param.T_ref)),
+            values.evaluate(param.D_p_dimensional(param.c_p_init(1), param.T_ref)),
             1 * 10 ** (-13),
             2,
         )