#933 remove reactions dict

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)
 

pybamm/models/full_battery_models/base_battery_model.py

@@ -196,8 +196,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
                     )

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/lithium_ion/base_lithium_ion_model.py

@@ -36,14 +36,10 @@ def set_standard_output_variables(self):
             }
         )
 
-    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
-        icd = " interfacial current density"
-        self.reactions = {
-            "main": {
-                "Negative": {"s": 1, "aj": "Negative electrode" + icd},
-                "Positive": {"s": 1, "aj": "Positive electrode" + icd},
-            }
-        }
+    def set_other_reaction_submodels_to_zero(self):
+        self.submodels["negative oxygen interface"] = pybamm.interface.NoReaction(
+            self.param, "Negative", "lithium-ion oxygen"
+        )
+        self.submodels["positive oxygen interface"] = pybamm.interface.NoReaction(
+            self.param, "Positive", "lithium-ion oxygen"
+        )

pybamm/models/full_battery_models/lithium_ion/dfn.py

@@ -38,6 +38,7 @@ def __init__(self, options=None, name="Doyle-Fuller-Newman model", build=True):
         self.set_tortuosity_submodels()
         self.set_convection_submodel()
         self.set_interfacial_submodel()
+        self.set_other_reaction_submodels_to_zero()
         self.set_particle_submodel()
         self.set_solid_submodel()
         self.set_electrolyte_submodel()

pybamm/models/full_battery_models/lithium_ion/spm.py

@@ -37,6 +37,7 @@ def __init__(self, options=None, name="Single Particle Model", build=True):
         self.set_tortuosity_submodels()
         self.set_convection_submodel()
         self.set_interfacial_submodel()
+        self.set_other_reaction_submodels_to_zero()
         self.set_particle_submodel()
         self.set_negative_electrode_submodel()
         self.set_electrolyte_submodel()

pybamm/models/full_battery_models/lithium_ion/spme.py

@@ -40,6 +40,7 @@ def __init__(
         self.set_tortuosity_submodels()
         self.set_convection_submodel()
         self.set_interfacial_submodel()
+        self.set_other_reaction_submodels_to_zero()
         self.set_particle_submodel()
         self.set_negative_electrode_submodel()
         self.set_electrolyte_submodel()

pybamm/models/submodels/base_submodel.py

@@ -47,12 +47,7 @@ class BaseSubModel:
     """
 
     def __init__(
-        self,
-        param,
-        domain=None,
-        reactions=None,
-        name="Unnamed submodel",
-        external=False,
+        self, param, domain=None, name="Unnamed submodel", external=False,
     ):
         super().__init__()
         self.param = param
@@ -67,7 +62,6 @@ def __init__(
 
         self.domain = domain
         self.set_domain_for_broadcast()
-        self.reactions = reactions
         self.name = name
 
         self.external = external

pybamm/models/submodels/electrode/base_electrode.py

@@ -19,8 +19,8 @@ class BaseElectrode(pybamm.BaseSubModel):
     **Extends:** :class:`pybamm.BaseSubModel`
     """
 
-    def __init__(self, param, domain, reactions=None, set_positive_potential=True):
-        super().__init__(param, domain, reactions)
+    def __init__(self, param, domain, set_positive_potential=True):
+        super().__init__(param, domain)
         self.set_positive_potential = set_positive_potential
 
     def _get_standard_potential_variables(self, phi_s):
@@ -190,4 +190,3 @@ def _get_standard_whole_cell_variables(self, variables):
             )
 
         return variables
-

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

@@ -20,8 +20,8 @@ class BaseModel(BaseElectrode):
     **Extends:** :class:`pybamm.electrode.BaseElectrode`
     """
 
-    def __init__(self, param, domain, reactions=None, set_positive_potential=True):
-        super().__init__(param, domain, reactions, set_positive_potential)
+    def __init__(self, param, domain, set_positive_potential=True):
+        super().__init__(param, domain, set_positive_potential)
 
     def set_boundary_conditions(self, variables):
 

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

@@ -19,8 +19,8 @@ class Full(BaseModel):
     **Extends:** :class:`pybamm.electrode.ohm.BaseModel`
     """
 
-    def __init__(self, param, domain, reactions):
-        super().__init__(param, domain, reactions)
+    def __init__(self, param, domain):
+        super().__init__(param, domain)
 
     def get_fundamental_variables(self):
 
@@ -59,10 +59,12 @@ def set_algebraic(self, variables):
 
         phi_s = variables[self.domain + " electrode potential"]
         i_s = variables[self.domain + " electrode current density"]
-        sum_j = sum(
-            variables[reaction[self.domain]["aj"]]
-            for reaction in self.reactions.values()
-        )
+
+        # All possible reactions. Some of these could be zero
+        j = variables[self.domain + " electrode interfacial current density"]
+        j_ox = variables[self.domain + " electrode oxygen interfacial current density"]
+
+        sum_j = j + j_ox
 
         self.algebraic[phi_s] = pybamm.div(i_s) + sum_j
 

pybamm/models/submodels/electrolyte_conductivity/base_electrolyte_conductivity.py

@@ -20,9 +20,8 @@ class BaseElectrolyteConductivity(pybamm.BaseSubModel):
     **Extends:** :class:`pybamm.BaseSubModel`
     """
 
-    def __init__(self, param, domain=None, reactions=None):
+    def __init__(self, param, domain=None):
         super().__init__(param, domain)
-        self.reactions = reactions
 
     def _get_standard_potential_variables(self, phi_e_n, phi_e_s, phi_e_p):
         """

pybamm/models/submodels/electrolyte_conductivity/full_conductivity.py

@@ -21,8 +21,8 @@ class Full(BaseElectrolyteConductivity):
     **Extends:** :class:`pybamm.electrolyte_conductivity.BaseElectrolyteConductivity`
     """
 
-    def __init__(self, param, reactions):
-        super().__init__(param, reactions=reactions)
+    def __init__(self, param):
+        super().__init__(param)
 
     def get_fundamental_variables(self):
         phi_e_n = pybamm.standard_variables.phi_e_n
@@ -51,14 +51,12 @@ def get_coupled_variables(self, variables):
     def set_algebraic(self, variables):
         phi_e = variables["Electrolyte potential"]
         i_e = variables["Electrolyte current density"]
-        sum_j = sum(
-            pybamm.Concatenation(
-                variables[reaction["Negative"]["aj"]],
-                pybamm.FullBroadcast(0, "separator", "current collector"),
-                variables[reaction["Positive"]["aj"]],
-            )
-            for reaction in self.reactions.values()
-        )
+
+        # All possible reactions. Some of these could be zero
+        j = variables["Interfacial current density"]
+        j_ox = variables["Oxygen interfacial current density"]
+
+        sum_j = j + j_ox
 
         self.algebraic = {phi_e: pybamm.div(i_e) - sum_j}
 

pybamm/models/submodels/electrolyte_conductivity/leading_order_conductivity.py

@@ -22,8 +22,8 @@ class LeadingOrder(BaseElectrolyteConductivity):
     **Extends:** :class:`pybamm.electrolyte_conductivity.BaseElectrolyteConductivity`
     """
 
-    def __init__(self, param, domain=None, reactions=None):
-        super().__init__(param, domain, reactions)
+    def __init__(self, param, domain=None):
+        super().__init__(param, domain)
 
     def get_coupled_variables(self, variables):
         ocp_n_av = variables["X-averaged negative electrode open circuit potential"]

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

@@ -22,8 +22,8 @@ class BaseModel(BaseElectrolyteConductivity):
     **Extends:** :class:`pybamm.electrolyte_conductivity.BaseElectrolyteConductivity`
     """
 
-    def __init__(self, param, domain, reactions):
-        super().__init__(param, domain, reactions)
+    def __init__(self, param, domain):
+        super().__init__(param, domain)
 
     def get_fundamental_variables(self):
         if self.domain == "Negative":
@@ -222,19 +222,21 @@ class FullAlgebraic(BaseModel):
      **Extends:** :class:`pybamm.electrolyte_conductivity.surface_potential_form.BaseFull`
     """  # noqa: E501
 
-    def __init__(self, param, domain, reactions):
-        super().__init__(param, domain, reactions)
+    def __init__(self, param, domain):
+        super().__init__(param, domain)
 
     def set_algebraic(self, variables):
         if self.domain == "Separator":
             return
 
         delta_phi = variables[self.domain + " electrode surface potential difference"]
         i_e = variables[self.domain + " electrolyte current density"]
-        sum_j = sum(
-            variables[reaction[self.domain]["aj"]]
-            for reaction in self.reactions.values()
-        )
+
+        # All possible reactions. Some of these could be zero
+        j = variables[self.domain + " electrode interfacial current density"]
+        j_ox = variables[self.domain + " electrode oxygen interfacial current density"]
+
+        sum_j = j + j_ox
         self.algebraic[delta_phi] = pybamm.div(i_e) - sum_j
 
 
@@ -253,8 +255,8 @@ class FullDifferential(BaseModel):
 
     """  # noqa: E501
 
-    def __init__(self, param, domain, reactions):
-        super().__init__(param, domain, reactions)
+    def __init__(self, param, domain):
+        super().__init__(param, domain)
 
     def set_rhs(self, variables):
         if self.domain == "Separator":
@@ -267,9 +269,11 @@ def set_rhs(self, variables):
 
         delta_phi = variables[self.domain + " electrode surface potential difference"]
         i_e = variables[self.domain + " electrolyte current density"]
-        sum_j = sum(
-            variables[reaction[self.domain]["aj"]]
-            for reaction in self.reactions.values()
-        )
+
+        # All possible reactions. Some of these could be zero
+        j = variables[self.domain + " electrode interfacial current density"]
+        j_ox = variables[self.domain + " electrode oxygen interfacial current density"]
+
+        sum_j = j + j_ox
 
         self.rhs[delta_phi] = 1 / C_dl * (pybamm.div(i_e) - sum_j)

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

@@ -23,8 +23,8 @@ class BaseLeadingOrderSurfaceForm(LeadingOrder):
     **Extends:** :class:`pybamm.electrolyte_conductivity.LeadingOrder`
     """
 
-    def __init__(self, param, domain, reactions):
-        super().__init__(param, domain, reactions)
+    def __init__(self, param, domain):
+        super().__init__(param, domain)
 
     def get_fundamental_variables(self):
 
@@ -94,19 +94,28 @@ class LeadingOrderDifferential(BaseLeadingOrderSurfaceForm):
 
     """
 
-    def __init__(self, param, domain, reactions):
-        super().__init__(param, domain, reactions)
+    def __init__(self, param, domain):
+        super().__init__(param, domain)
 
     def set_rhs(self, variables):
         if self.domain == "Separator":
             return
 
         param = self.param
 
-        sum_j = sum(
-            variables["X-averaged " + reaction[self.domain]["aj"].lower()]
-            for reaction in self.reactions.values()
-        )
+        # All possible reactions. Some of these could be zero
+        j = variables[
+            "X-averaged "
+            + self.domain.lower()
+            + " electrode interfacial current density"
+        ]
+        j_ox = variables[
+            "X-averaged "
+            + self.domain.lower()
+            + " electrode oxygen interfacial current density"
+        ]
+
+        sum_j = j + j_ox
 
         sum_j_av = variables[
             "X-averaged "
@@ -141,17 +150,26 @@ class LeadingOrderAlgebraic(BaseLeadingOrderSurfaceForm):
     **Extends:** :class:`BaseLeadingOrderSurfaceForm`
     """
 
-    def __init__(self, param, domain, reactions):
-        super().__init__(param, domain, reactions)
+    def __init__(self, param, domain):
+        super().__init__(param, domain)
 
     def set_algebraic(self, variables):
         if self.domain == "Separator":
             return
 
-        sum_j = sum(
-            variables["X-averaged " + reaction[self.domain]["aj"].lower()]
-            for reaction in self.reactions.values()
-        )
+        # All possible reactions. Some of these could be zero
+        j = variables[
+            "X-averaged "
+            + self.domain.lower()
+            + " electrode interfacial current density"
+        ]
+        j_ox = variables[
+            "X-averaged "
+            + self.domain.lower()
+            + " electrode oxygen interfacial current density"
+        ]
+
+        sum_j = j + j_ox
 
         sum_j_av = variables[
             "X-averaged "

pybamm/models/submodels/electrolyte_diffusion/base_electrolyte_diffusion.py

@@ -17,8 +17,8 @@ class BaseElectrolyteDiffusion(pybamm.BaseSubModel):
     **Extends:** :class:`pybamm.BaseSubModel`
     """
 
-    def __init__(self, param, reactions=None):
-        super().__init__(param, reactions=reactions)
+    def __init__(self, param):
+        super().__init__(param)
 
     def _get_standard_concentration_variables(self, c_e_n, c_e_s, c_e_p):
         """