#492 reorganise and add power

Author: valentinsulzer

examples/scripts/compare_lithium_ion.py

@@ -16,7 +16,7 @@
     pybamm.set_logging_level("INFO")
 
 # load models
-options = {"thermal": "isothermal", "operating mode": "voltage"}
+options = {"thermal": "isothermal", "operating mode": "power"}
 models = [
     pybamm.lithium_ion.SPM(options),
     pybamm.lithium_ion.SPMe(options),
@@ -27,7 +27,7 @@
 # load parameter values and process models and geometry
 param = models[0].default_parameter_values
 param["Typical current [A]"] = 1.0
-param["Voltage function"] = 4  # voltage
+param["Power function"] = -1  # voltage
 for model in models:
     param.process_model(model)
 
@@ -51,5 +51,15 @@
     solutions[i] = model.default_solver.solve(model, t_eval)
 
 # plot
-plot = pybamm.QuickPlot(models, mesh, solutions)
+output_variables = [
+    "Negative particle surface concentration",
+    "Electrolyte concentration",
+    "Positive particle surface concentration",
+    "Current [A]",
+    "Negative electrode potential [V]",
+    "Electrolyte potential [V]",
+    "Terminal power [W]",
+    "Terminal voltage [V]",
+]
+plot = pybamm.QuickPlot(models, mesh, solutions, output_variables)
 plot.dynamic_plot()

pybamm/models/base_model.py

@@ -388,13 +388,21 @@ def check_algebraic_equations(self, post_discretisation):
         After discretisation, there must be at least one StateVector in each algebraic
         equation
         """
+        vars_in_bcs = set()
+        for var, side_eqn in self.boundary_conditions.items():
+            for side, (eqn, typ) in side_eqn.items():
+                vars_in_bcs.update(
+                    [x.id for x in eqn.pre_order() if isinstance(x, pybamm.Variable)]
+                )
         if not post_discretisation:
             # After the model has been defined, each algebraic equation key should
-            # appear in that algebraic equation
+            # appear in that algebraic equation, or in the boundary conditions
             # this has been relaxed for concatenations for now
             for var, eqn in self.algebraic.items():
-                if not any(x.id == var.id for x in eqn.pre_order()) and not isinstance(
-                    var, pybamm.Concatenation
+                if not (
+                    any(x.id == var.id for x in eqn.pre_order())
+                    or var.id in vars_in_bcs
+                    or isinstance(var, pybamm.Concatenation)
                 ):
                     raise pybamm.ModelError(
                         "each variable in the algebraic eqn keys must appear in the eqn"

pybamm/models/full_battery_models/base_battery_model.py

@@ -162,7 +162,12 @@ def options(self, extra_options):
                     raise pybamm.OptionError("option {} not recognised".format(name))
 
         # Some standard checks to make sure options are compatible
-        if options["operating mode"] not in ["current", "voltage"]:
+        if options["operating mode"] not in [
+            "current",
+            "voltage",
+            "power",
+            "arbitrary",
+        ]:
             raise pybamm.OptionError(
                 "operating mode '{}' not recognised".format(options["operating mode"])
             )
@@ -491,14 +496,22 @@ def set_external_circuit_submodel(self):
         e.g. (not necessarily constant-) current, voltage, etc
         """
         if self.options["operating mode"] == "current":
+            self.submodels["external circuit"] = pybamm.external_circuit.CurrentControl(
+                self.param
+            )
+        elif self.options["operating mode"] == "voltage":
+            self.submodels["external circuit"] = pybamm.external_circuit.VoltageControl(
+                self.param
+            )
+        elif self.options["operating mode"] == "power":
+            self.submodels["external circuit"] = pybamm.external_circuit.PowerControl(
+                self.param
+            )
+        elif self.options["operating mode"] == "arbitrary":
             self.submodels[
                 "external circuit"
-            ] = pybamm.external_circuit.CurrentControl(self.param)
-        if self.options["operating mode"] == "voltage":
-            self.submodels[
-                "external circuit"
-            ] = pybamm.external_circuit.VoltageControl(self.param)
-            
+            ] = pybamm.external_circuit.FunctionControl(self.param)
+
     def set_tortuosity_submodels(self):
         self.submodels["electrolyte tortuosity"] = pybamm.tortuosity.Bruggeman(
             self.param, "Electrolyte"
@@ -732,6 +745,10 @@ def set_voltage_variables(self):
         self.events["Minimum voltage"] = voltage - self.param.voltage_low_cut
         self.events["Maximum voltage"] = voltage - self.param.voltage_high_cut
 
+        # Power
+        I_dim = self.variables["Current [A]"]
+        self.variables.update({"Terminal power [W]": I_dim * V_dim})
+
     def set_soc_variables(self):
         """
         Set variables relating to the state of charge.

pybamm/models/full_battery_models/lead_acid/base_lead_acid_model.py

@@ -53,19 +53,6 @@ def default_solver(self):
 
     def set_standard_output_variables(self):
         super().set_standard_output_variables()
-        # Current
-        i_cell = pybamm.standard_parameters_lead_acid.current_with_time
-        i_cell_dim = (
-            pybamm.standard_parameters_lead_acid.dimensional_current_density_with_time
-        )
-        I = pybamm.standard_parameters_lead_acid.dimensional_current_with_time
-        self.variables.update(
-            {
-                "Total current density": i_cell,
-                "Total current density [A.m-2]": i_cell_dim,
-                "Current [A]": I,
-            }
-        )
 
         # Time
         time_scale = pybamm.standard_parameters_lead_acid.tau_discharge

pybamm/models/submodels/external_circuit/__init__.py

@@ -1,3 +1,6 @@
 from .base_external_circuit import BaseModel
 from .current_control_external_circuit import CurrentControl
+from .function_control_external_circuit import FunctionControl
 from .voltage_control_external_circuit import VoltageControl
+from .power_control_external_circuit import PowerControl
+

pybamm/models/submodels/external_circuit/function_control_external_circuit.py

@@ -0,0 +1,48 @@
+#
+# External circuit with an arbitrary function
+#
+import pybamm
+from .base_external_circuit import BaseModel
+
+
+class FunctionControl(BaseModel):
+    """External circuit with an arbitrary function. """
+
+    def __init__(self, param):
+        super().__init__(param)
+
+    def external_circuit_function(self, I, V):
+        """
+        Function that fixes the current, I (in Amps), or voltage, V (in Volts), or a
+        combination of the two
+        """
+        return pybamm.FunctionParameter("External circuit function", I, V)
+
+    def get_fundamental_variables(self):
+        # Current is a variable
+        param = self.param
+        i_cell = pybamm.Variable("Total current density")
+        I = i_cell * abs(param.I_typ)
+        i_cell_dim = I / (param.n_electrodes_parallel * param.A_cc)
+
+        variables = {
+            "Total current density": i_cell,
+            "Total current density [A.m-2]": i_cell_dim,
+            "Current [A]": I,
+        }
+
+        return variables
+
+    def set_initial_conditions(self, variables):
+        # Initial condition as a guess for consistent initial conditions
+        i_cell = variables["Total current density"]
+        self.initial_conditions[i_cell] = self.param.current_with_time
+
+    def set_algebraic(self, variables):
+        # External circuit submodels are always equations on the current
+        # The external circuit function should fix either the current, or the voltage,
+        # or a combination (e.g. I*V for power control)
+        i_cell = variables["Total current density"]
+        I = variables["Current [A]"]
+        V = variables["Terminal voltage [V]"]
+        self.algebraic[i_cell] = self.external_circuit_function(I, V)

pybamm/models/submodels/external_circuit/power_control_external_circuit.py

@@ -0,0 +1,18 @@
+#
+# External circuit with power control
+#
+import pybamm
+from .function_control_external_circuit import FunctionControl
+
+
+class PowerControl(FunctionControl):
+    """External circuit with power control. """
+
+    def __init__(self, param):
+        super().__init__(param)
+
+    def external_circuit_function(self, I, V):
+        return I * V - pybamm.FunctionParameter(
+            "Power function", pybamm.t * self.param.timescale
+        )
+

pybamm/models/submodels/external_circuit/voltage_control_external_circuit.py

@@ -2,39 +2,17 @@
 # External circuit with voltage control
 #
 import pybamm
-from .base_external_circuit import BaseModel
+from .function_control_external_circuit import FunctionControl
 
 
-class VoltageControl(BaseModel):
+class VoltageControl(FunctionControl):
     """External circuit with voltage control. """
 
     def __init__(self, param):
         super().__init__(param)
 
-    def get_fundamental_variables(self):
-        # Current is a variable
-        param = self.param
-        i_cell = pybamm.Variable("Total current density")
-        I = i_cell * abs(param.I_typ)
-        i_cell_dim = I / (param.n_electrodes_parallel * param.A_cc)
+    def external_circuit_function(self, I, V):
+        return V - pybamm.FunctionParameter(
+            "Voltage function", pybamm.t * self.param.timescale
+        )
 
-        variables = {
-            "Total current density": i_cell,
-            "Total current density [A.m-2]": i_cell_dim,
-            "Current [A]": I,
-        }
-
-        return variables
-
-    def set_initial_conditions(self, variables):
-        # Initial condition as a guess for consistent initial conditions
-        i_cell = variables["Total current density"]
-        self.initial_conditions[i_cell] = self.param.current_with_time
-
-    def set_algebraic(self, variables):
-        # External circuit submodels are always equations on the current
-        # Fix voltage to be equal to terminal voltage
-        i_cell = variables["Total current density"]
-        V = variables["Terminal voltage"]
-        # TODO: find a way to get rid of 0 * i_cell
-        self.algebraic[i_cell] = V - self.param.voltage_with_time + 0 * i_cell

pybamm/parameters/standard_parameters_lead_acid.py

@@ -495,11 +495,3 @@ def U_p(c_e_p, T):
     dimensional_current_with_time / I_typ * pybamm.Function(np.sign, I_typ)
 )
 
-# Voltage
-dimensional_voltage_with_time = pybamm.FunctionParameter(
-    "Voltage function", pybamm.t * timescale
-)
-voltage_with_time = (
-    dimensional_voltage_with_time - (U_p_ref - U_n_ref)
-) / potential_scale
-

pybamm/parameters/standard_parameters_lithium_ion.py

@@ -452,10 +452,3 @@ def dUdT_p(c_s_p):
     dimensional_current_with_time / I_typ * pybamm.Function(np.sign, I_typ)
 )
 
-# Voltage
-dimensional_voltage_with_time = pybamm.FunctionParameter(
-    "Voltage function", pybamm.t * timescale
-)
-voltage_with_time = (
-    dimensional_voltage_with_time - (U_p_ref - U_n_ref)
-) / potential_scale