#492 get potentiostatic working for 1D lithium-ion

Author: valentinsulzer

examples/scripts/compare_lead_acid.py

@@ -16,16 +16,23 @@
     pybamm.set_logging_level("INFO")
 
 # load models
+options = {"operating mode": "voltage"}
 models = [
-    pybamm.lead_acid.LOQS(),
-    pybamm.lead_acid.FOQS(),
-    pybamm.lead_acid.Composite(),
-    pybamm.lead_acid.Full(),
+    pybamm.lead_acid.LOQS(options),
+    # pybamm.lead_acid.FOQS(options),
+    # pybamm.lead_acid.Composite(options),
+    pybamm.lead_acid.Full(options),
 ]
 
 # load parameter values and process models and geometry
 param = models[0].default_parameter_values
-param.update({"Typical current [A]": 10, "Initial State of Charge": 1})
+param.update(
+    {
+        "Typical current [A]": 10,
+        "Initial State of Charge": 1,
+        "Voltage function": 14 / 6,
+    }
+)
 for model in models:
     param.process_model(model)
 

examples/scripts/compare_lithium_ion.py

@@ -18,17 +18,16 @@
 # load models
 options = {"thermal": "isothermal", "operating mode": "voltage"}
 models = [
-    # pybamm.lithium_ion.SPM(options),
-    # pybamm.lithium_ion.SPMe(options),
-    # pybamm.lithium_ion.DFN(),
-    pybamm.lithium_ion.DFN(options)
+    pybamm.lithium_ion.SPM(options),
+    pybamm.lithium_ion.SPMe(options),
+    pybamm.lithium_ion.DFN(options),
 ]
 
 
 # load parameter values and process models and geometry
 param = models[0].default_parameter_values
 param["Typical current [A]"] = 1.0
-param["Voltage function"] = 4
+param["Voltage function"] = 4  # voltage
 for model in models:
     param.process_model(model)
 

pybamm/models/base_model.py

@@ -338,7 +338,8 @@ def check_well_determined(self, post_discretisation):
         vars_in_rhs_keys = set()
         vars_in_algebraic_keys = set()
         vars_in_eqns = set()
-        # Get all variables ids from rhs and algebraic keys and equations
+        # Get all variables ids from rhs and algebraic keys and equations, and
+        # from boundary conditions
         # For equations we look through the whole expression tree.
         # "Variables" can be Concatenations so we also have to look in the whole
         # expression tree
@@ -356,11 +357,16 @@ def check_well_determined(self, post_discretisation):
             vars_in_eqns.update(
                 [x.id for x in eqn.pre_order() if isinstance(x, pybamm.Variable)]
             )
+        for var, side_eqn in self.boundary_conditions.items():
+            for side, (eqn, typ) in side_eqn.items():
+                vars_in_eqns.update(
+                    [x.id for x in eqn.pre_order() if isinstance(x, pybamm.Variable)]
+                )
         # If any keys are repeated between rhs and algebraic then the model is
         # overdetermined
         if not set(vars_in_rhs_keys).isdisjoint(vars_in_algebraic_keys):
             raise pybamm.ModelError("model is overdetermined (repeated keys)")
-        # If any algebraic keys don't appear in the eqns then the model is
+        # If any algebraic keys don't appear in the eqns (or bcs) then the model is
         # overdetermined (but rhs keys can be absent from the eqns, e.g. dcdt = -1 is
         # fine)
         # Skip this step after discretisation, as any variables in the equations will

pybamm/models/full_battery_models/base_battery_model.py

@@ -466,10 +466,10 @@ def build_model(self):
             )
             self.update(submodel)
 
-        pybamm.logger.debug("Setting voltage variables")
+        pybamm.logger.debug("Setting voltage variables ({})".format(self.name))
         self.set_voltage_variables()
 
-        pybamm.logger.debug("Setting SoC variables")
+        pybamm.logger.debug("Setting SoC variables ({})".format(self.name))
         self.set_soc_variables()
 
         # Massive hack for consistent delta_phi = phi_s - phi_e with SPMe

pybamm/models/full_battery_models/lithium_ion/dfn.py

@@ -85,10 +85,10 @@ def set_particle_submodel(self):
     def set_solid_submodel(self):
 
         self.submodels["negative electrode"] = pybamm.electrode.ohm.Full(
-            self.param, "Negative", self.reactions, self.options["operating mode"]
+            self.param, "Negative", self.reactions
         )
         self.submodels["positive electrode"] = pybamm.electrode.ohm.Full(
-            self.param, "Positive", self.reactions, self.options["operating mode"]
+            self.param, "Positive", self.reactions
         )
 
     def set_electrolyte_submodel(self):

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

@@ -19,9 +19,8 @@ class Full(BaseModel):
     **Extends:** :class:`pybamm.electrode.ohm.BaseModel`
     """
 
-    def __init__(self, param, domain, reactions, operating_mode):
+    def __init__(self, param, domain, reactions):
         super().__init__(param, domain, reactions)
-        self.operating_mode = operating_mode
 
     def get_fundamental_variables(self):
 
@@ -80,14 +79,11 @@ def set_boundary_conditions(self, variables):
 
         elif self.domain == "Positive":
             lbc = (pybamm.Scalar(0), "Neumann")
-            if self.operating_mode == "current":
-                sigma_eff = self.param.sigma_p * tor
-                rbc = (
-                    i_boundary_cc / pybamm.boundary_value(-sigma_eff, "right"),
-                    "Neumann",
-                )
-            elif self.operating_mode == "voltage":
-                rbc = (self.param.voltage_with_time, "Dirichlet")
+            sigma_eff = self.param.sigma_p * tor
+            rbc = (
+                i_boundary_cc / pybamm.boundary_value(-sigma_eff, "right"),
+                "Neumann",
+            )
 
         self.boundary_conditions[phi_s] = {"left": lbc, "right": rbc}
 

pybamm/models/submodels/external_circuit/voltage_control_external_circuit.py

@@ -32,17 +32,10 @@ def set_initial_conditions(self, variables):
         self.initial_conditions[i_cell] = self.param.current_with_time
 
     def set_algebraic(self, variables):
-        # Read off the current from the solid current density
-        # TODO: investigate whether defining the current differently gives better
-        # results (e.g. based on electrolyte current density in the separator)
-        # This also needs to be defined in a model-inpedendent way (currently assumes
-        # Ohm's law)
+        # External circuit submodels are always equations on the current
+        # Fix voltage to be equal to terminal voltage
         i_cell = variables["Total current density"]
-        phi_s_p = variables["Positive electrode potential"]
-        tor_p = variables["Positive electrode tortuosity"]
-        sigma_eff = self.param.sigma_p * tor_p
-        i_s_p_right = -pybamm.boundary_value(
-            sigma_eff, "right"
-        ) * pybamm.BoundaryGradient(phi_s_p, "right")
-
-        self.algebraic[i_cell] = i_cell - i_s_p_right
+        # TODO: change this to terminal voltage
+        V = variables["Local current collector potential difference"]
+        # TODO: find a way to get rid of 0 * i_cell
+        self.algebraic[i_cell] = V - self.param.voltage_with_time + 0 * i_cell