Merge pull request pybamm-team#2552 from pybamm-team/issue-2551-electrolyte

Issue 2551 electrolyte

Author: valentinsulzer

CHANGELOG.md

@@ -9,6 +9,7 @@
 
 ## Bug fixes
 
+- Fixed electrolyte conservation when options {"surface form": "algebraic"} are used
 - Fixed "constant concentration" electrolyte model so that "porosity times concentration" is conserved when porosity changes ([#2529](https://github.com/pybamm-team/PyBaMM/pull/2529))
 - Fix installation on `Google Colab` (`pybtex` and `Colab` issue) ([#2526](https://github.com/pybamm-team/PyBaMM/pull/2526))
 

pybamm/expression_tree/unary_operators.py

@@ -1252,6 +1252,14 @@ def boundary_value(symbol, side):
         return BoundaryValue(symbol, side)
 
 
+def boundary_gradient(symbol, side):
+    # Gradient of a broadcast is zero
+    if isinstance(symbol, pybamm.Broadcast):
+        return 0 * symbol.reduce_one_dimension()
+    else:
+        return BoundaryGradient(symbol, side)
+
+
 def sign(symbol):
     """Returns a :class:`Sign` object."""
     if isinstance(symbol, pybamm.Broadcast):

pybamm/models/submodels/electrolyte_conductivity/surface_potential_form/full_surface_form_conductivity.py

@@ -60,6 +60,9 @@ def get_coupled_variables(self, variables):
                 + i_boundary_cc / sigma_eff
             )
             variables[f"{Domain} electrolyte current density"] = i_e
+            variables[
+                f"Divergence of {self.domain} electrolyte current density"
+            ] = pybamm.div(i_e)
 
             phi_s = variables[f"{Domain} electrode potential"]
             phi_e = phi_s - delta_phi
@@ -91,8 +94,8 @@ def get_coupled_variables(self, variables):
 
             # Update boundary conditions (for indefinite integral)
             self.boundary_conditions[c_e_s] = {
-                "left": (pybamm.BoundaryGradient(c_e_s, "left"), "Neumann"),
-                "right": (pybamm.BoundaryGradient(c_e_s, "right"), "Neumann"),
+                "left": (pybamm.boundary_gradient(c_e_s, "left"), "Neumann"),
+                "right": (pybamm.boundary_gradient(c_e_s, "right"), "Neumann"),
             }
 
         variables[f"{Domain} electrolyte potential"] = phi_e
@@ -111,6 +114,34 @@ def get_coupled_variables(self, variables):
             variables.update(self._get_standard_current_variables(i_e))
             variables.update(self._get_electrolyte_overpotentials(variables))
 
+        # save boundary conditons as variables
+        if self.domain == "negative":
+            grad_c_e = pybamm.boundary_gradient(c_e, "right")
+            grad_left = -i_boundary_cc * pybamm.boundary_value(1 / sigma_eff, "left")
+            grad_right = (
+                (i_boundary_cc / pybamm.boundary_value(conductivity, "right"))
+                - pybamm.boundary_value(param.chiRT_over_Fc(c_e, T), "right") * grad_c_e
+                - i_boundary_cc * pybamm.boundary_value(1 / sigma_eff, "right")
+            )
+
+        elif self.domain == "positive":
+            T = variables["Positive electrode temperature"]
+            grad_c_e = pybamm.boundary_gradient(c_e, "left")
+            grad_left = (
+                (i_boundary_cc / pybamm.boundary_value(conductivity, "left"))
+                - pybamm.boundary_value(param.chiRT_over_Fc(c_e, T), "left") * grad_c_e
+                - i_boundary_cc * pybamm.boundary_value(1 / sigma_eff, "left")
+            )
+            grad_right = -i_boundary_cc * pybamm.boundary_value(1 / sigma_eff, "right")
+
+        if self.domain in ["negative", "positive"]:
+            variables.update(
+                {
+                    f"{self.domain} grad(delta_phi) left": grad_left,
+                    f"{self.domain} grad(delta_phi) right": grad_right,
+                    f"{self.domain} grad(c_e) internal": grad_c_e,
+                }
+            )
         return variables
 
     def _get_conductivities(self, variables):
@@ -146,41 +177,20 @@ def set_boundary_conditions(self, variables):
         if self.domain == "separator":
             return
 
-        param = self.param
-
-        conductivity, sigma_eff = self._get_conductivities(variables)
-        i_boundary_cc = variables["Current collector current density"]
         c_e = variables[f"{Domain} electrolyte concentration"]
         delta_phi = variables[f"{Domain} electrode surface potential difference"]
 
-        if self.domain == "negative":
-            T = variables["Negative electrode temperature"]
-            c_e_flux = pybamm.BoundaryGradient(c_e, "right")
-            flux_left = -i_boundary_cc * pybamm.BoundaryValue(1 / sigma_eff, "left")
-            flux_right = (
-                (i_boundary_cc / pybamm.BoundaryValue(conductivity, "right"))
-                - pybamm.BoundaryValue(param.chiRT_over_Fc(c_e, T), "right") * c_e_flux
-                - i_boundary_cc * pybamm.BoundaryValue(1 / sigma_eff, "right")
-            )
+        grad_left = variables[f"{self.domain} grad(delta_phi) left"]
+        grad_right = variables[f"{self.domain} grad(delta_phi) right"]
+        grad_c_e = variables[f"{self.domain} grad(c_e) internal"]
 
-            lbc = (flux_left, "Neumann")
-            rbc = (flux_right, "Neumann")
+        lbc = (grad_left, "Neumann")
+        rbc = (grad_right, "Neumann")
+        if self.domain == "negative":
             lbc_c_e = (pybamm.Scalar(0), "Neumann")
-            rbc_c_e = (c_e_flux, "Neumann")
-
+            rbc_c_e = (grad_c_e, "Neumann")
         elif self.domain == "positive":
-            T = variables["Positive electrode temperature"]
-            c_e_flux = pybamm.BoundaryGradient(c_e, "left")
-            flux_left = (
-                (i_boundary_cc / pybamm.BoundaryValue(conductivity, "left"))
-                - pybamm.BoundaryValue(param.chiRT_over_Fc(c_e, T), "left") * c_e_flux
-                - i_boundary_cc * pybamm.BoundaryValue(1 / sigma_eff, "left")
-            )
-            flux_right = -i_boundary_cc * pybamm.BoundaryValue(1 / sigma_eff, "right")
-
-            lbc = (flux_left, "Neumann")
-            rbc = (flux_right, "Neumann")
-            lbc_c_e = (c_e_flux, "Neumann")
+            lbc_c_e = (grad_c_e, "Neumann")
             rbc_c_e = (pybamm.Scalar(0), "Neumann")
 
         # TODO: check why we still need the boundary conditions for c_e, once we have

pybamm/models/submodels/electrolyte_diffusion/full_diffusion.py

@@ -76,6 +76,13 @@ def get_coupled_variables(self, variables):
         N_e = N_e_diffusion + N_e_migration + N_e_convection
 
         variables.update(self._get_standard_flux_variables(N_e))
+        variables.update(
+            {
+                "Electrolyte diffusion flux": N_e_diffusion,
+                "Electrolyte migration flux": N_e_migration,
+                "Electrolyte convection flux": N_e_convection,
+            }
+        )
 
         return variables
 
@@ -85,12 +92,16 @@ def set_rhs(self, variables):
 
         eps_c_e = variables["Porosity times concentration"]
         c_e = variables["Electrolyte concentration"]
-        N_e = variables["Electrolyte flux"]
+        T = variables["Cell temperature"]
+        N_e_diffusion = variables["Electrolyte diffusion flux"]
+        N_e_convection = variables["Electrolyte convection flux"]
+        N_e = N_e_diffusion + N_e_convection
         div_Vbox = variables["Transverse volume-averaged acceleration"]
 
         sum_s_j = variables["Sum of electrolyte reaction source terms"]
+        sum_a_j = variables["Sum of volumetric interfacial current densities"]
         sum_s_j.print_name = "a"
-        source_terms = sum_s_j / self.param.gamma_e
+        source_terms = (sum_s_j - param.t_plus(c_e, T) * sum_a_j) / self.param.gamma_e
 
         self.rhs = {
             eps_c_e: -pybamm.div(N_e) / param.C_e + source_terms - c_e * div_Vbox

pybamm/models/submodels/interface/kinetics/diffusion_limited.py

@@ -84,7 +84,7 @@ def _get_diffusion_limited_current_density(self, variables):
                 N_ox_neg_sep_interface = (
                     -pybamm.boundary_value(tor_s, "left")
                     * param.curlyD_ox
-                    * pybamm.BoundaryGradient(c_ox_s, "left")
+                    * pybamm.boundary_gradient(c_ox_s, "left")
                 )
                 N_ox_neg_sep_interface.domains = {"primary": "current collector"}
 

pybamm/models/submodels/thermal/pouch_cell/pouch_cell_2D_current_collectors.py

@@ -115,8 +115,8 @@ def set_boundary_conditions(self, variables):
             / self.param.delta
         )
 
-        T_av_n = pybamm.BoundaryValue(T_av, "negative tab")
-        T_av_p = pybamm.BoundaryValue(T_av, "positive tab")
+        T_av_n = pybamm.boundary_value(T_av, "negative tab")
+        T_av_p = pybamm.boundary_value(T_av, "positive tab")
 
         self.boundary_conditions = {
             T_av: {

pybamm/spatial_methods/finite_volume.py

@@ -615,8 +615,6 @@ def add_ghost_nodes(self, symbol, discretised_symbol, bcs):
         n = submesh.npts
         second_dim_repeats = self._get_auxiliary_domain_repeats(symbol.domains)
 
-        bcs_vector = pybamm.Vector([])  # starts empty
-
         lbc_value, lbc_type = bcs["left"]
         rbc_value, rbc_type = bcs["right"]
 
@@ -748,23 +746,23 @@ def add_neumann_values(self, symbol, discretised_gradient, bcs, domain):
             n_bcs += 1
 
         # Add any values from Neumann boundary conditions to the bcs vector
-        if lbc_type == "Neumann":
+        if lbc_type == "Neumann" and lbc_value != 0:
             lbc_sub_matrix = coo_matrix(([1], ([0], [0])), shape=(n + n_bcs, 1))
             lbc_matrix = csr_matrix(kron(eye(second_dim_repeats), lbc_sub_matrix))
             if lbc_value.evaluates_to_number():
                 left_bc = lbc_value * pybamm.Vector(np.ones(second_dim_repeats))
             else:
                 left_bc = lbc_value
             lbc_vector = pybamm.Matrix(lbc_matrix) @ left_bc
-        elif lbc_type == "Dirichlet":
+        elif lbc_type == "Dirichlet" or (lbc_type == "Neumann" and lbc_value == 0):
             lbc_vector = pybamm.Vector(np.zeros((n + n_bcs) * second_dim_repeats))
         else:
             raise ValueError(
                 "boundary condition must be Dirichlet or Neumann, not '{}'".format(
-                    lbc_type
+                    rbc_type
                 )
             )
-        if rbc_type == "Neumann":
+        if rbc_type == "Neumann" and rbc_value != 0:
             rbc_sub_matrix = coo_matrix(
                 ([1], ([n + n_bcs - 1], [0])), shape=(n + n_bcs, 1)
             )
@@ -774,7 +772,7 @@ def add_neumann_values(self, symbol, discretised_gradient, bcs, domain):
             else:
                 right_bc = rbc_value
             rbc_vector = pybamm.Matrix(rbc_matrix) @ right_bc
-        elif rbc_type == "Dirichlet":
+        elif rbc_type == "Dirichlet" or (rbc_type == "Neumann" and rbc_value == 0):
             rbc_vector = pybamm.Vector(np.zeros((n + n_bcs) * second_dim_repeats))
         else:
             raise ValueError(

tests/integration/test_models/standard_output_tests.py

@@ -508,7 +508,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])
-        np.testing.assert_array_almost_equal(diff, 0)
+        if self.model.options["surface form"] == "differential":
+            np.testing.assert_array_almost_equal(diff, 0, decimal=9)
+        else:
+            np.testing.assert_array_almost_equal(diff, 0, decimal=14)
 
     def test_concentration_profile(self):
         """Test continuity of the concentration profile. Test average concentration is
@@ -560,11 +563,14 @@ def test_splitting(self):
 
     def test_all(self):
         self.test_concentration_limit()
-        self.test_conservation()
         self.test_concentration_profile()
         self.test_fluxes()
         self.test_splitting()
 
+        if isinstance(self.model, pybamm.lithium_ion.BaseModel):
+            # electrolyte is not conserved in lead-acid models
+            self.test_conservation()
+
 
 class PotentialTests(BaseOutputTest):
     def __init__(self, model, param, disc, solution, operating_condition):

tests/unit/test_expression_tree/test_unary_operators.py

@@ -577,6 +577,15 @@ def test_boundary_value(self):
         ):
             pybamm.boundary_value(symbol_on_edges, "right")
 
+    def test_boundary_gradient(self):
+        var = pybamm.Variable("var", domain=["negative electrode"])
+        grad = pybamm.boundary_gradient(var, "right")
+        self.assertIsInstance(grad, pybamm.BoundaryGradient)
+
+        zero = pybamm.PrimaryBroadcast(0, ["negative electrode"])
+        grad = pybamm.boundary_gradient(zero, "right")
+        self.assertEqual(grad, 0)
+
     def test_unary_simplifications(self):
         a = pybamm.Scalar(0)
         b = pybamm.Scalar(1)