Issue 1903 power

CHANGELOG.md

@@ -3,6 +3,7 @@
 
 ## Features
 
+-   Added new models for power control and resistance control ([#1917](https://github.com/pybamm-team/PyBaMM/pull/1917))
 -   Initial concentrations can now be provided as a function of `r` as well as `x` ([#1866](https://github.com/pybamm-team/PyBaMM/pull/1866))
 
 ## Breaking changes

docs/source/models/submodels/external_circuit/explicit_control_external_circuit.rst

@@ -0,0 +1,13 @@
+Explicit control external circuit
+=================================
+
+Current is explicitly specified, either by a function or in terms of other variables.
+
+.. autoclass:: pybamm.external_circuit.ExplicitCurrentControl
+    :members:
+
+.. autoclass:: pybamm.external_circuit.ExplicitPowerControl
+    :members:
+
+.. autoclass:: pybamm.external_circuit.ExplicitResistanceControl
+    :members:

docs/source/models/submodels/external_circuit/function_control_external_circuit.rst

@@ -10,5 +10,8 @@ Function control external circuit
 .. autoclass:: pybamm.external_circuit.PowerFunctionControl
     :members:
 
+.. autoclass:: pybamm.external_circuit.ResistanceFunctionControl
+    :members:
+
 .. autoclass:: pybamm.external_circuit.CCCVFunctionControl
     :members:

docs/source/models/submodels/external_circuit/index.rst

@@ -11,5 +11,5 @@ variable to be constant.
 .. toctree::
   :maxdepth: 1
 
-  current_control_external_circuit
+  explicit_control_external_circuit
   function_control_external_circuit

examples/notebooks/models/using-submodels.ipynb

@@ -335,7 +335,7 @@
    "source": [
     "Submodels can be added to the `model.submodels` dictionary in the same way that we changed the submodels earlier. \n",
     "\n",
-    "We use the simplest model for the external circuit, which is the \"current control\" submodel"
+    "We use the simplest model for the external circuit, which is the explicit \"current control\" submodel"
    ]
   },
   {
@@ -344,7 +344,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "model.submodels[\"external circuit\"] = pybamm.external_circuit.CurrentControl(model.param)"
+    "model.submodels[\"external circuit\"] = pybamm.external_circuit.ExplicitCurrentControl(model.param)"
    ]
   },
   {

examples/scripts/DFN.py

@@ -5,46 +5,41 @@
 import pybamm
 import numpy as np
 
-pybamm.set_logging_level("INFO")
-
-# load model
-model = pybamm.lithium_ion.DFN()
-# create geometry
-geometry = model.default_geometry
-
-# load parameter values and process model and geometry
-param = model.default_parameter_values
-param.process_geometry(geometry)
-param.process_model(model)
-
-# set mesh
-var = pybamm.standard_spatial_vars
-var_pts = {var.x_n: 30, var.x_s: 30, var.x_p: 30, var.r_n: 10, var.r_p: 10}
-mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts)
-
-# discretise model
-disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
-disc.process_model(model)
+models = [
+    pybamm.lithium_ion.DFN()  # {"operating mode": "explicit power"}),
+    # pybamm.lithium_ion.DFN({"operating mode": "power"}),
+    # pybamm.lithium_ion.DFN({"operating mode": "differential power"}),
+]
+
+# set parameters and discretise models
+for i, model in enumerate(models):
+    # create geometry
+    params = model.default_parameter_values
+    # params.update({"Power function [W]": 4}, check_already_exists=False)
+    params.update({"Current function [A]": 0})
+    geometry = model.default_geometry
+    params.process_model(model)
+    params.process_geometry(geometry)
+    mesh = pybamm.Mesh(geometry, model.default_submesh_types, model.default_var_pts)
+    disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
+    disc.process_model(model)
 
 # solve model
+# pybamm.set_logging_level("DEBUG")
+solutions = [None] * len(models)
 t_eval = np.linspace(0, 3600, 100)
-solver = pybamm.CasadiSolver(mode="safe", atol=1e-6, rtol=1e-3)
-solution = solver.solve(model, t_eval)
-
-# plot
-plot = pybamm.QuickPlot(
-    solution,
+for i, model in enumerate(models):
+    print(model.name)
+    # solver = pybamm.ScipySolver()
+    solver = pybamm.CasadiSolver()
+    solutions[i] = solver.solve(model, t_eval)
+    # print(solutions[i].solve_time)
+    # print(solutions[i].integration_time)
+pybamm.dynamic_plot(
+    solutions,
     [
-        "Negative particle concentration [mol.m-3]",
-        "Electrolyte concentration [mol.m-3]",
-        "Positive particle concentration [mol.m-3]",
-        "Current [A]",
-        "Negative electrode potential [V]",
-        "Electrolyte potential [V]",
-        "Positive electrode potential [V]",
         "Terminal voltage [V]",
+        "Current [A]",
+        "Terminal power [W]",
     ],
-    time_unit="seconds",
-    spatial_unit="um",
 )
-plot.dynamic_plot()

examples/scripts/custom_model.py

@@ -11,7 +11,7 @@
 model = pybamm.lithium_ion.BaseModel(name="my li-ion model")
 
 # set choice of submodels
-model.submodels["external circuit"] = pybamm.external_circuit.CurrentControl(
+model.submodels["external circuit"] = pybamm.external_circuit.ExplicitCurrentControl(
     model.param
 )
 model.submodels["current collector"] = pybamm.current_collector.Uniform(model.param)

pybamm/models/full_battery_models/base_battery_model.py

@@ -65,12 +65,19 @@ class BatteryModelOptions(pybamm.FuzzyDict):
                 A 2-tuple can be provided for different behaviour in negative and
                 positive electrodes.
             * "operating mode" : str
-                Sets the operating mode for the model. Can be "current" (default),
-                "voltage" or "power". Alternatively, the operating mode can be
-                controlled with an arbitrary function by passing the function directly
-                as the option. In this case the function must define the residual of
-                an algebraic equation. The applied current will be solved for such
-                that the algebraic constraint is satisfied.
+                Sets the operating mode for the model. This determines how the current
+                is set. Can be:
+
+                - "current" (default) : the current is explicity supplied
+                - "voltage"/"power"/"resistance" : solve an algebraic equation for \
+                    current such that voltage/power/resistance is correct
+                - "differential power"/"differential resistance" : solve a \
+                    differential equation for the power or resistance
+                - "explicit power"/"explicit resistance" : current is defined in terms \
+                    of the voltage such that power/resistance is correct
+                - callable : if a callable is given as this option, the function \
+                    defines the residual of an algebraic equation. The applied current \
+                    will be solved for such that the algebraic constraint is satisfied.
             * "particle" : str
                 Sets the submodel to use to describe behaviour within the particle.
                 Can be "Fickian diffusion" (default), "uniform profile",
@@ -186,7 +193,17 @@ def __init__(self, extra_options):
                 "reaction-driven",
                 "stress and reaction-driven",
             ],
-            "operating mode": ["current", "voltage", "power", "CCCV"],
+            "operating mode": [
+                "current",
+                "voltage",
+                "power",
+                "differential power",
+                "explicit power",
+                "resistance",
+                "differential resistance",
+                "explicit resistance",
+                "CCCV",
+            ],
             "particle": [
                 "Fickian diffusion",
                 "fast diffusion",
@@ -885,27 +902,36 @@ 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
-            )
+            model = pybamm.external_circuit.ExplicitCurrentControl(self.param)
         elif self.options["operating mode"] == "voltage":
-            self.submodels[
-                "external circuit"
-            ] = pybamm.external_circuit.VoltageFunctionControl(self.param)
+            model = pybamm.external_circuit.VoltageFunctionControl(self.param)
         elif self.options["operating mode"] == "power":
-            self.submodels[
-                "external circuit"
-            ] = pybamm.external_circuit.PowerFunctionControl(self.param)
+            model = pybamm.external_circuit.PowerFunctionControl(
+                self.param, "algebraic"
+            )
+        elif self.options["operating mode"] == "differential power":
+            model = pybamm.external_circuit.PowerFunctionControl(
+                self.param, "differential without max"
+            )
+        elif self.options["operating mode"] == "explicit power":
+            model = pybamm.external_circuit.ExplicitPowerControl(self.param)
+        elif self.options["operating mode"] == "resistance":
+            model = pybamm.external_circuit.ResistanceFunctionControl(
+                self.param, "algebraic"
+            )
+        elif self.options["operating mode"] == "differential resistance":
+            model = pybamm.external_circuit.ResistanceFunctionControl(
+                self.param, "differential without max"
+            )
+        elif self.options["operating mode"] == "explicit resistance":
+            model = pybamm.external_circuit.ExplicitResistanceControl(self.param)
         elif self.options["operating mode"] == "CCCV":
-            self.submodels[
-                "external circuit"
-            ] = pybamm.external_circuit.CCCVFunctionControl(self.param)
+            model = pybamm.external_circuit.CCCVFunctionControl(self.param)
         elif callable(self.options["operating mode"]):
-            self.submodels[
-                "external circuit"
-            ] = pybamm.external_circuit.FunctionControl(
+            model = pybamm.external_circuit.FunctionControl(
                 self.param, self.options["operating mode"]
             )
+        self.submodels["external circuit"] = model
 
     def set_tortuosity_submodels(self):
         self.submodels["electrolyte tortuosity"] = pybamm.tortuosity.Bruggeman(
@@ -1125,10 +1151,11 @@ def set_voltage_variables(self):
         # Hack to avoid division by zero if i_cc is exactly zero
         # If i_cc is zero, i_cc_not_zero becomes 1. But multiplying by sign(i_cc) makes
         # the local resistance 'zero' (really, it's not defined when i_cc is zero)
-        i_cc_not_zero = ((i_cc > 0) + (i_cc < 0)) * i_cc + (i_cc >= 0) * (i_cc <= 0)
-        i_cc_dim_not_zero = ((i_cc_dim > 0) + (i_cc_dim < 0)) * i_cc_dim + (
-            i_cc_dim >= 0
-        ) * (i_cc_dim <= 0)
+        def x_not_zero(x):
+            return ((x > 0) + (x < 0)) * x + (x >= 0) * (x <= 0)
+
+        i_cc_not_zero = x_not_zero(i_cc)
+        i_cc_dim_not_zero = x_not_zero(i_cc_dim)
 
         self.variables.update(
             {
@@ -1179,9 +1206,16 @@ def set_voltage_variables(self):
             )
         )
 
-        # Power
+        # Power and resistance
         I_dim = self.variables["Current [A]"]
-        self.variables.update({"Terminal power [W]": I_dim * V_dim})
+        I_dim_not_zero = x_not_zero(I_dim)
+        self.variables.update(
+            {
+                "Terminal power [W]": I_dim * V_dim,
+                "Power [W]": I_dim * V_dim,
+                "Resistance [Ohm]": pybamm.sign(I_dim) * V_dim / I_dim_not_zero,
+            }
+        )
 
     def set_degradation_variables(self):
         """

pybamm/models/full_battery_models/lead_acid/loqs.py

@@ -64,7 +64,7 @@ def set_external_circuit_submodel(self):
         if self.options["operating mode"] == "current":
             self.submodels[
                 "leading order external circuit"
-            ] = pybamm.external_circuit.LeadingOrderCurrentControl(self.param)
+            ] = pybamm.external_circuit.LeadingOrderExplicitCurrentControl(self.param)
         elif self.options["operating mode"] == "voltage":
             self.submodels[
                 "leading order external circuit"

pybamm/models/submodels/current_collector/homogeneous_current_collector.py

@@ -21,17 +21,20 @@ class Uniform(BaseModel):
     def __init__(self, param):
         super().__init__(param)
 
+    def get_fundamental_variables(self):
+        phi_s_cn = pybamm.PrimaryBroadcast(0, "current collector")
+        variables = self._get_standard_negative_potential_variables(phi_s_cn)
+        return variables
+
     def get_coupled_variables(self, variables):
 
         # TODO: grad not implemented for 2D yet
         i_cc = pybamm.Scalar(0)
         i_boundary_cc = pybamm.PrimaryBroadcast(
             variables["Total current density"], "current collector"
         )
-        phi_s_cn = pybamm.PrimaryBroadcast(0, "current collector")
 
-        variables = self._get_standard_negative_potential_variables(phi_s_cn)
-        variables.update(self._get_standard_current_variables(i_cc, i_boundary_cc))
+        variables = self._get_standard_current_variables(i_cc, i_boundary_cc)
 
         # Hack to get the leading-order current collector current density
         # Note that this should be different from the actual (composite) current

pybamm/models/submodels/electrolyte_diffusion/full_diffusion.py

@@ -113,13 +113,15 @@ def set_boundary_conditions(self, variables):
             T = variables["Cell temperature"]
             tor = variables["Electrolyte tortuosity"]
             i_boundary_cc = variables["Current collector current density"]
-            dce_dx = (
-                -(1 - param.t_plus(c_e, T))
+            lbc = (
+                pybamm.boundary_value(
+                    -(1 - param.t_plus(c_e, T))
+                    / (tor * param.gamma_e * param.D_e(c_e, T)),
+                    "left",
+                )
                 * i_boundary_cc
                 * param.C_e
-                / (tor * param.gamma_e * param.D_e(c_e, T))
             )
-            lbc = pybamm.boundary_value(dce_dx, "left")
         else:
             # left bc at anode/current collector interface
             lbc = pybamm.Scalar(0)

pybamm/models/submodels/external_circuit/__init__.py

@@ -1,9 +1,15 @@
 from .base_external_circuit import BaseModel, LeadingOrderBaseModel
-from .current_control_external_circuit import CurrentControl, LeadingOrderCurrentControl
+from .explicit_control_external_circuit import (
+    ExplicitCurrentControl,
+    ExplicitPowerControl,
+    ExplicitResistanceControl,
+    LeadingOrderExplicitCurrentControl,
+)
 from .function_control_external_circuit import (
     FunctionControl,
     VoltageFunctionControl,
     PowerFunctionControl,
+    ResistanceFunctionControl,
     CCCVFunctionControl,
     LeadingOrderFunctionControl,
     LeadingOrderVoltageFunctionControl,