Merge pull request #948 from pybamm-team/issue-933-remove-reactions-dict

Issue 933 remove reactions dict

CHANGELOG.md

@@ -11,6 +11,7 @@
 
 ## Optimizations
 
+-   Added some extra simplifications to the expression tree ([#971](https://github.com/pybamm-team/PyBaMM/pull/971))
 -   Changed the behaviour of "safe" mode in `CasadiSolver` ([#956](https://github.com/pybamm-team/PyBaMM/pull/956))
 -   Sped up model building ([#927](https://github.com/pybamm-team/PyBaMM/pull/927))
 -   Changed default solver for lead-acid to `CasadiSolver` ([#927](https://github.com/pybamm-team/PyBaMM/pull/927))
@@ -29,6 +30,7 @@
 
 ## Breaking changes
 
+-   Changed the implementation of reactions in submodels ([#948](https://github.com/pybamm-team/PyBaMM/pull/948))
 -   Removed some inputs like `T_inf`, `R_g` and activation energies to some of the standard function parameters. This is because each of those inputs is specific to a particular function (e.g. the reference temperature at which the function was measured). To change a property such as the activation energy, users should create a new function, specifying the relevant property as a `Parameter` or `InputParameter` ([#942](https://github.com/pybamm-team/PyBaMM/pull/942))
 -   The thermal option 'xyz-lumped' has been removed. The option 'thermal current collector' has also been removed ([#938](https://github.com/pybamm-team/PyBaMM/pull/938))
 -   The 'C-rate' parameter has been deprecated. Use 'Current function [A]' instead. The cell capacity can be accessed as 'Cell capacity [A.h]', and used to calculate current from C-rate ([#952](https://github.com/pybamm-team/PyBaMM/pull/952))

examples/scripts/compare_lead_acid.py

@@ -25,7 +25,7 @@
 
 # load parameter values and process models and geometry
 param = models[0].default_parameter_values
-param.update({"Current function [A]": 85, "Initial State of Charge": 1})
+param.update({"Current function [A]": 17, "Initial State of Charge": 1})
 for model in models:
     param.process_model(model)
 

examples/scripts/compare_lithium_ion_particle_distribution.py

@@ -61,7 +61,7 @@ def positive_distribution(x):
 solutions = [None] * len(models)
 t_eval = np.linspace(0, 3600, 100)
 for i, model in enumerate(models):
-    solutions[i] = model.default_solver.solve(model, t_eval)
+    solutions[i] = pybamm.CasadiSolver().solve(model, t_eval)
 
 
 output_variables = [

pybamm/input/parameters/lead-acid/anodes/lead_Sulzer2019/parameters.csv

@@ -18,18 +18,18 @@ Negative electrode morphological parameter,0.6,srinivasan2003mathematical,
 Negative electrode capacity [C.m-3],3473000000,,
 ,,,
 # Interfacial reactions,,,
-Negative electrode cation signed stoichiometry,-1,,
+Negative electrode cation signed stoichiometry,1,,
 Negative electrode electrons in reaction,2,,
 Negative electrode reference exchange-current density [A.m-2],0.06,srinivasan2003mathematical,
-Signed stoichiometry of cations (oxygen reaction),-4,,
-Signed stoichiometry of water (oxygen reaction),1,,
-Signed stoichiometry of oxygen (oxygen reaction),-1,,
+Signed stoichiometry of cations (oxygen reaction),4,,
+Signed stoichiometry of water (oxygen reaction),-1,,
+Signed stoichiometry of oxygen (oxygen reaction),1,,
 Electrons in oxygen reaction,4,,
 Negative electrode reference exchange-current density (oxygen) [A.m-2],2.5E-32,srinivasan2003mathematical,
 Reference oxygen molecule concentration [mol.m-3],1000,srinivasan2003mathematical,
 Oxygen reference OCP vs SHE [V],1.229,srinivasan2003mathematical,
-Signed stoichiometry of cations (hydrogen reaction),-2,,
-Signed stoichiometry of hydrogen (hydrogen reaction),1,,
+Signed stoichiometry of cations (hydrogen reaction),2,,
+Signed stoichiometry of hydrogen (hydrogen reaction),-1,,
 Electrons in hydrogen reaction,2,,
 Negative electrode reference exchange-current density (hydrogen) [A.m-2],1.56E-11,srinivasan2003mathematical,
 Hydrogen reference OCP vs SHE [V],0,srinivasan2003mathematical,

pybamm/input/parameters/lead-acid/cathodes/lead_dioxide_Sulzer2019/parameters.csv

@@ -18,18 +18,18 @@ Positive electrode morphological parameter,0.6,srinivasan2003mathematical,
 Positive electrode capacity [C.m-3],2745000000,,
 ,,,
 # Interfacial reactions,,,
-Positive electrode cation signed stoichiometry,-3,,
+Positive electrode cation signed stoichiometry,3,,
 Positive electrode electrons in reaction,2,,
 Positive electrode reference exchange-current density [A.m-2],0.004,srinivasan2003mathematical,
-Signed stoichiometry of cations (oxygen reaction),-4,,
-Signed stoichiometry of water (oxygen reaction),1,,
-Signed stoichiometry of oxygen (oxygen reaction),-1,,
+Signed stoichiometry of cations (oxygen reaction),4,,
+Signed stoichiometry of water (oxygen reaction),-1,,
+Signed stoichiometry of oxygen (oxygen reaction),1,,
 Electrons in oxygen reaction,4,,
 Positive electrode reference exchange-current density (oxygen) [A.m-2],2.5E-23,srinivasan2003mathematical,
 Reference oxygen molecule concentration [mol.m-3],1000,srinivasan2003mathematical,
 Oxygen reference OCP vs SHE [V],1.229,srinivasan2003mathematical,
-Signed stoichiometry of cations (hydrogen reaction),-2,,
-Signed stoichiometry of hydrogen (hydrogen reaction),1,,
+Signed stoichiometry of cations (hydrogen reaction),2,,
+Signed stoichiometry of hydrogen (hydrogen reaction),-1,,
 Electrons in hydrogen reaction,2,,
 Positive electrode reference exchange-current density (hydrogen) [A.m-2],0,srinivasan2003mathematical,
 Hydrogen reference OCP vs SHE [V],0,srinivasan2003mathematical,

pybamm/models/full_battery_models/base_battery_model.py

@@ -186,8 +186,7 @@ def options(self, extra_options):
         ):
             if len(options["side reactions"]) > 0:
                 raise pybamm.OptionError(
-                    """
-                    must use surface formulation to solve {!s} with side reactions
+                    """must use surface formulation to solve {!s} with side reactions
                     """.format(
                         self
                     )
@@ -283,6 +282,24 @@ def set_standard_output_variables(self):
                 {"y": var.y, "y [m]": var.y * L_y, "z": var.z, "z [m]": var.z * L_z}
             )
 
+        # Initialize "total reaction" variables
+        self.variables.update(
+            {
+                "Sum of electrolyte reaction source terms": 0,
+                "Sum of negative electrode electrolyte reaction source terms": 0,
+                "Sum of positive electrode electrolyte reaction source terms": 0,
+                "Sum of x-averaged negative electrode "
+                "electrolyte reaction source terms": 0,
+                "Sum of x-averaged positive electrode "
+                "electrolyte reaction source terms": 0,
+                "Sum of interfacial current densities": 0,
+                "Sum of negative electrode interfacial current densities": 0,
+                "Sum of positive electrode interfacial current densities": 0,
+                "Sum of x-averaged negative electrode interfacial current densities": 0,
+                "Sum of x-averaged positive electrode interfacial current densities": 0,
+            }
+        )
+
     def build_fundamental_and_external(self):
         # Get the fundamental variables
         for submodel_name, submodel in self.submodels.items():
@@ -339,11 +356,11 @@ def build_coupled_variables(self):
                         if len(submodels) == 1 or count == 100:
                             # no more submodels to try
                             raise pybamm.ModelError(
-                                """Submodel "{}" requires the variable {}, but it cannot be found.
-                                Check the selected submodels provide all of the required
-                                variables.""".format(
+                                "Submodel '{}' requires the variable {}, ".format(
                                     submodel_name, key
                                 )
+                                + "but it cannot be found. Check the selected "
+                                "submodels provide all of the required variables."
                             )
                         else:
                             # try setting coupled variables on next loop through

pybamm/models/full_battery_models/lead_acid/base_lead_acid_model.py

@@ -53,34 +53,6 @@ def default_solver(self):
         else:
             return pybamm.CasadiSolver(mode="safe")
 
-    def set_reactions(self):
-
-        # Should probably refactor as this is a bit clunky at the moment
-        # Maybe each reaction as a Reaction class so we can just list names of classes
-        param = self.param
-        icd = " interfacial current density"
-        self.reactions = {
-            "main": {
-                "Negative": {"s": -param.s_plus_n_S, "aj": "Negative electrode" + icd},
-                "Positive": {"s": -param.s_plus_p_S, "aj": "Positive electrode" + icd},
-            }
-        }
-        if "oxygen" in self.options["side reactions"]:
-            self.reactions["oxygen"] = {
-                "Negative": {
-                    "s": -param.s_plus_Ox,
-                    "s_ox": -param.s_ox_Ox,
-                    "aj": "Negative electrode oxygen" + icd,
-                },
-                "Positive": {
-                    "s": -param.s_plus_Ox,
-                    "s_ox": -param.s_ox_Ox,
-                    "aj": "Positive electrode oxygen" + icd,
-                },
-            }
-            self.reactions["main"]["Negative"]["s_ox"] = 0
-            self.reactions["main"]["Positive"]["s_ox"] = 0
-
     def set_soc_variables(self):
         "Set variables relating to the state of charge."
         # State of Charge defined as function of dimensionless electrolyte concentration

pybamm/models/full_battery_models/lead_acid/basic_full.py

@@ -28,7 +28,7 @@ class BasicFull(BaseModel):
     **Extends:** :class:`pybamm.lead_acid.BaseModel`
     """
 
-    def __init__(self, name="Full model"):
+    def __init__(self, name="Basic full model"):
         super().__init__({}, name)
         # `param` is a class containing all the relevant parameters and functions for
         # this model. These are purely symbolic at this stage, and will be set by the
@@ -244,15 +244,19 @@ def __init__(self, name="Full model"):
         ######################
         # Electrolyte concentration
         ######################
-        N_e = -tor * param.D_e(c_e, T) * pybamm.grad(c_e) + param.C_e * c_e * v
+        N_e = (
+            -tor * param.D_e(c_e, T) * pybamm.grad(c_e)
+            + param.C_e * param.t_plus(c_e) * i_e / param.gamma_e
+            + param.C_e * c_e * v
+        )
         s = pybamm.Concatenation(
-            -pybamm.PrimaryBroadcast(param.s_plus_n_S, "negative electrode"),
+            pybamm.PrimaryBroadcast(param.s_plus_n_S, "negative electrode"),
             pybamm.PrimaryBroadcast(0, "separator"),
-            -pybamm.PrimaryBroadcast(param.s_plus_p_S, "positive electrode"),
+            pybamm.PrimaryBroadcast(param.s_plus_p_S, "positive electrode"),
         )
         self.rhs[c_e] = (1 / eps) * (
             -pybamm.div(N_e) / param.C_e
-            + (s - param.t_plus(c_e)) * j / param.gamma_e
+            + s * j / param.gamma_e
             - c_e * deps_dt
             - c_e * div_V
         )
@@ -286,6 +290,7 @@ def __init__(self, name="Full model"):
             "Terminal voltage [V]": param.U_p_ref - param.U_n_ref + pot * voltage,
             "x [m]": pybamm.standard_spatial_vars.x * param.L_x,
             "x": pybamm.standard_spatial_vars.x,
+            "Porosity": eps,
             "Volume-averaged velocity": v,
             "X-averaged separator transverse volume-averaged velocity": div_V_s,
         }

pybamm/models/full_battery_models/lead_acid/full.py

@@ -35,7 +35,6 @@ def __init__(self, options=None, name="Full model", build=True):
         super().__init__(options, name)
 
         self.set_external_circuit_submodel()
-        self.set_reactions()
         self.set_interfacial_submodel()
         self.set_porosity_submodel()
         self.set_tortuosity_submodels()
@@ -85,8 +84,8 @@ def set_interfacial_submodel(self):
 
     def set_solid_submodel(self):
         if self.options["surface form"] is False:
-            submod_n = pybamm.electrode.ohm.Full(self.param, "Negative", self.reactions)
-            submod_p = pybamm.electrode.ohm.Full(self.param, "Positive", self.reactions)
+            submod_n = pybamm.electrode.ohm.Full(self.param, "Negative")
+            submod_p = pybamm.electrode.ohm.Full(self.param, "Positive")
         else:
             submod_n = pybamm.electrode.ohm.SurfaceForm(self.param, "Negative")
             submod_p = pybamm.electrode.ohm.SurfaceForm(self.param, "Positive")
@@ -99,28 +98,28 @@ def set_electrolyte_submodel(self):
         surf_form = pybamm.electrolyte_conductivity.surface_potential_form
 
         self.submodels["electrolyte diffusion"] = pybamm.electrolyte_diffusion.Full(
-            self.param, self.reactions
+            self.param
         )
 
         if self.options["surface form"] is False:
             self.submodels[
                 "electrolyte conductivity"
-            ] = pybamm.electrolyte_conductivity.Full(self.param, self.reactions)
+            ] = pybamm.electrolyte_conductivity.Full(self.param)
         elif self.options["surface form"] == "differential":
             for domain in ["Negative", "Separator", "Positive"]:
                 self.submodels[
                     domain.lower() + " electrolyte conductivity"
-                ] = surf_form.FullDifferential(self.param, domain, self.reactions)
+                ] = surf_form.FullDifferential(self.param, domain)
         elif self.options["surface form"] == "algebraic":
             for domain in ["Negative", "Separator", "Positive"]:
                 self.submodels[
                     domain.lower() + " electrolyte conductivity"
-                ] = surf_form.FullAlgebraic(self.param, domain, self.reactions)
+                ] = surf_form.FullAlgebraic(self.param, domain)
 
     def set_side_reaction_submodels(self):
         if "oxygen" in self.options["side reactions"]:
             self.submodels["oxygen diffusion"] = pybamm.oxygen_diffusion.Full(
-                self.param, self.reactions
+                self.param
             )
             self.submodels["positive oxygen interface"] = pybamm.interface.ForwardTafel(
                 self.param, "Positive", "lead-acid oxygen"

pybamm/models/full_battery_models/lead_acid/higher_order.py

@@ -36,7 +36,6 @@ def __init__(self, options=None, name="Composite model", build=True):
 
         self.set_external_circuit_submodel()
         self.set_leading_order_model()
-        self.set_reactions()
         # Electrolyte submodel to get first-order concentrations
         self.set_electrolyte_diffusion_submodel()
         self.set_other_species_diffusion_submodels()
@@ -98,6 +97,24 @@ def set_leading_order_model(self):
             leading_order_model.variables["X-averaged electrolyte concentration"]
         ]
 
+        # Reset sums
+        self.variables.update(
+            {
+                "Sum of electrolyte reaction source terms": 0,
+                "Sum of negative electrode electrolyte reaction source terms": 0,
+                "Sum of positive electrode electrolyte reaction source terms": 0,
+                "Sum of x-averaged negative electrode "
+                "electrolyte reaction source terms": 0,
+                "Sum of x-averaged positive electrode "
+                "electrolyte reaction source terms": 0,
+                "Sum of interfacial current densities": 0,
+                "Sum of negative electrode interfacial current densities": 0,
+                "Sum of positive electrode interfacial current densities": 0,
+                "Sum of x-averaged negative electrode interfacial current densities": 0,
+                "Sum of x-averaged positive electrode interfacial current densities": 0,
+            }
+        )
+
     def set_average_interfacial_submodel(self):
         self.submodels[
             "x-averaged negative interface"
@@ -204,12 +221,12 @@ def __init__(self, options=None, name="FOQS model", build=True):
     def set_electrolyte_diffusion_submodel(self):
         self.submodels[
             "electrolyte diffusion"
-        ] = pybamm.electrolyte_diffusion.FirstOrder(self.param, self.reactions)
+        ] = pybamm.electrolyte_diffusion.FirstOrder(self.param)
 
     def set_other_species_diffusion_submodels(self):
         if "oxygen" in self.options["side reactions"]:
             self.submodels["oxygen diffusion"] = pybamm.oxygen_diffusion.FirstOrder(
-                self.param, self.reactions
+                self.param
             )
 
     def set_full_porosity_submodel(self):
@@ -234,12 +251,12 @@ def __init__(self, options=None, name="Composite model", build=True):
     def set_electrolyte_diffusion_submodel(self):
         self.submodels[
             "electrolyte diffusion"
-        ] = pybamm.electrolyte_diffusion.Composite(self.param, self.reactions)
+        ] = pybamm.electrolyte_diffusion.Composite(self.param)
 
     def set_other_species_diffusion_submodels(self):
         if "oxygen" in self.options["side reactions"]:
             self.submodels["oxygen diffusion"] = pybamm.oxygen_diffusion.Composite(
-                self.param, self.reactions
+                self.param
             )
 
     def set_full_porosity_submodel(self):
@@ -278,14 +295,12 @@ def __init__(
     def set_electrolyte_diffusion_submodel(self):
         self.submodels[
             "electrolyte diffusion"
-        ] = pybamm.electrolyte_diffusion.Composite(
-            self.param, self.reactions, extended="distributed"
-        )
+        ] = pybamm.electrolyte_diffusion.Composite(self.param, extended="distributed")
 
     def set_other_species_diffusion_submodels(self):
         if "oxygen" in self.options["side reactions"]:
             self.submodels["oxygen diffusion"] = pybamm.oxygen_diffusion.Composite(
-                self.param, self.reactions, extended="distributed"
+                self.param, extended="distributed"
             )
 
 
@@ -302,12 +317,10 @@ def __init__(self, options=None, name="Extended composite model (average)"):
     def set_electrolyte_diffusion_submodel(self):
         self.submodels[
             "electrolyte diffusion"
-        ] = pybamm.electrolyte_diffusion.Composite(
-            self.param, self.reactions, extended="average"
-        )
+        ] = pybamm.electrolyte_diffusion.Composite(self.param, extended="average")
 
     def set_other_species_diffusion_submodels(self):
         if "oxygen" in self.options["side reactions"]:
             self.submodels["oxygen diffusion"] = pybamm.oxygen_diffusion.Composite(
-                self.param, self.reactions, extended="average"
+                self.param, extended="average"
             )

pybamm/models/full_battery_models/lead_acid/loqs.py

@@ -34,7 +34,6 @@ def __init__(self, options=None, name="LOQS model", build=True):
         super().__init__(options, name)
 
         self.set_external_circuit_submodel()
-        self.set_reactions()
         self.set_interfacial_submodel()
         self.set_convection_submodel()
         self.set_porosity_submodel()
@@ -185,25 +184,23 @@ def set_electrolyte_submodel(self):
             for domain in ["Negative", "Separator", "Positive"]:
                 self.submodels[
                     "leading-order " + domain.lower() + " electrolyte conductivity"
-                ] = surf_form.LeadingOrderDifferential(
-                    self.param, domain, self.reactions
-                )
+                ] = surf_form.LeadingOrderDifferential(self.param, domain)
 
         elif self.options["surface form"] == "algebraic":
             for domain in ["Negative", "Separator", "Positive"]:
                 self.submodels[
                     "leading-order " + domain.lower() + " electrolyte conductivity"
-                ] = surf_form.LeadingOrderAlgebraic(self.param, domain, self.reactions)
+                ] = surf_form.LeadingOrderAlgebraic(self.param, domain)
 
         self.submodels[
             "electrolyte diffusion"
-        ] = pybamm.electrolyte_diffusion.LeadingOrder(self.param, self.reactions)
+        ] = pybamm.electrolyte_diffusion.LeadingOrder(self.param)
 
     def set_side_reaction_submodels(self):
         if "oxygen" in self.options["side reactions"]:
             self.submodels[
                 "leading-order oxygen diffusion"
-            ] = pybamm.oxygen_diffusion.LeadingOrder(self.param, self.reactions)
+            ] = pybamm.oxygen_diffusion.LeadingOrder(self.param)
             self.submodels[
                 "leading-order positive oxygen interface"
             ] = pybamm.interface.ForwardTafel(