Repository: pybamm-team/PyBaMM

Issue:
Add Discharged energy [W.h] as variable
**Summary**
Add Discharged energy [W.h] as variable

 **Motivation**
Like discharged capacity but using current * voltage. Useful for performance metrics and minimal overhead


Discussion:

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 22e31f970..52450aac8 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,7 @@
 
 ## Features
 
+-   Added "Discharge energy [W.h]", which is the integral of the power in Watts, as an optional output. Set the option "calculate discharge energy" to "true" to get this output ("false" by default, since it can slow down some of the simple models) ([#1969](https://github.com/pybamm-team/PyBaMM/pull/1969)))
 -   Added an option "calculate heat source for isothermal models" to choose whether or not the heat generation terms are computed when running models with the option `thermal="isothermal"`  ([#1958](https://github.com/pybamm-team/PyBaMM/pull/1958))
 
 ## Bug fixes
diff --git a/examples/notebooks/models/using-submodels.ipynb b/examples/notebooks/models/using-submodels.ipynb
index 655801716..923a2675e 100644
--- a/examples/notebooks/models/using-submodels.ipynb
+++ b/examples/notebooks/models/using-submodels.ipynb
@@ -344,7 +344,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "model.submodels[\"external circuit\"] = pybamm.external_circuit.ExplicitCurrentControl(model.param)"
+    "model.submodels[\"external circuit\"] = pybamm.external_circuit.ExplicitCurrentControl(model.param, model.options)"
    ]
   },
   {
diff --git a/examples/scripts/custom_model.py b/examples/scripts/custom_model.py
index 6121b3fb5..0703d0c35 100644
--- a/examples/scripts/custom_model.py
+++ b/examples/scripts/custom_model.py
@@ -12,7 +12,7 @@ model = pybamm.lithium_ion.BaseModel(name="my li-ion model")
 
 # set choice of submodels
 model.submodels["external circuit"] = pybamm.external_circuit.ExplicitCurrentControl(
-    model.param
+    model.param, model.options
 )
 model.submodels["current collector"] = pybamm.current_collector.Uniform(model.param)
 model.submodels["thermal"] = pybamm.thermal.isothermal.Isothermal(model.param)
diff --git a/pybamm/models/full_battery_models/base_battery_model.py b/pybamm/models/full_battery_models/base_battery_model.py
index dc128b438..d9be896f2 100644
--- a/pybamm/models/full_battery_models/base_battery_model.py
+++ b/pybamm/models/full_battery_models/base_battery_model.py
@@ -20,6 +20,10 @@ class BatteryModelOptions(pybamm.FuzzyDict):
         where a 2-tuple of strings can be provided instead to indicate a different
         option for the negative and positive electrodes.
 
+            * "calculate discharge energy": str
+                Whether to calculate the discharge energy. Must be one of "true" or
+                "false". "false" is the default, since calculating the discharge
+                energy can be computationally expensive for simple models like SPM.
             * "cell geometry" : str
                 Sets the geometry of the cell. Can be "pouch" (default) or
                 "arbitrary". The arbitrary geometry option solves a 1D electrochemical
@@ -168,6 +172,7 @@ class BatteryModelOptions(pybamm.FuzzyDict):
 
     def __init__(self, extra_options):
         self.possible_options = {
+            "calculate discharge energy": ["false", "true"],
             "cell geometry": ["arbitrary", "pouch"],
             "calculate heat source for isothermal models": ["false", "true"],
             "convection": ["none", "uniform transverse", "full transverse"],
@@ -922,34 +927,44 @@ class BaseBatteryModel(pybamm.BaseModel):
         e.g. (not necessarily constant-) current, voltage, etc
         """
         if self.options["operating mode"] == "current":
-            model = pybamm.external_circuit.ExplicitCurrentControl(self.param)
+            model = pybamm.external_circuit.ExplicitCurrentControl(
+                self.param, self.options
+            )
         elif self.options["operating mode"] == "voltage":
-            model = pybamm.external_circuit.VoltageFunctionControl(self.param)
+            model = pybamm.external_circuit.VoltageFunctionControl(
+                self.param, self.options
+            )
         elif self.options["operating mode"] == "power":
             model = pybamm.external_circuit.PowerFunctionControl(
-                self.param, "algebraic"
+                self.param, self.options, "algebraic"
             )
         elif self.options["operating mode"] == "differential power":
             model = pybamm.external_circuit.PowerFunctionControl(
-                self.param, "differential without max"
+                self.param, self.options, "differential without max"
             )
         elif self.options["operating mode"] == "explicit power":
-            model = pybamm.external_circuit.ExplicitPowerControl(self.param)
+            model = pybamm.external_circuit.ExplicitPowerControl(
+                self.param, self.options
+            )
         elif self.options["operating mode"] == "resistance":
             model = pybamm.external_circuit.ResistanceFunctionControl(
-                self.param, "algebraic"
+                self.param, self.options, "algebraic"
             )
         elif self.options["operating mode"] == "differential resistance":
             model = pybamm.external_circuit.ResistanceFunctionControl(
-                self.param, "differential without max"
+                self.param, self.options, "differential without max"
             )
         elif self.options["operating mode"] == "explicit resistance":
-            model = pybamm.external_circuit.ExplicitResistanceControl(self.param)
+            model = pybamm.external_circuit.ExplicitResistanceControl(
+                self.param, self.options
+            )
         elif self.options["operating mode"] == "CCCV":
-            model = pybamm.external_circuit.CCCVFunctionControl(self.param)
+            model = pybamm.external_circuit.CCCVFunctionControl(
+                self.param, self.options
+            )
         elif callable(self.options["operating mode"]):
             model = pybamm.external_circuit.FunctionControl(
-                self.param, self.options["operating mode"]
+                self.param, self.options["operating mode"], self.options
             )
         self.submodels["external circuit"] = model
 
diff --git a/pybamm/models/full_battery_models/lead_acid/loqs.py b/pybamm/models/full_battery_models/lead_acid/loqs.py
index 602fe4d7f..d5ccd7c1d 100644
--- a/pybamm/models/full_battery_models/lead_acid/loqs.py
+++ b/pybamm/models/full_battery_models/lead_acid/loqs.py
@@ -64,20 +64,26 @@ class LOQS(BaseModel):
         if self.options["operating mode"] == "current":
             self.submodels[
                 "leading order external circuit"
-            ] = pybamm.external_circuit.LeadingOrderExplicitCurrentControl(self.param)
+            ] = pybamm.external_circuit.LeadingOrderExplicitCurrentControl(
+                self.param, self.options
+            )
         elif self.options["operating mode"] == "voltage":
             self.submodels[
                 "leading order external circuit"
-            ] = pybamm.external_circuit.LeadingOrderVoltageFunctionControl(self.param)
+            ] = pybamm.external_circuit.LeadingOrderVoltageFunctionControl(
+                self.param, self.options
+            )
         elif self.options["operating mode"] == "power":
             self.submodels[
                 "leading order external circuit"
-            ] = pybamm.external_circuit.LeadingOrderPowerFunctionControl(self.param)
+            ] = pybamm.external_circuit.LeadingOrderPowerFunctionControl(
+                self.param, self.options
+            )
         elif callable(self.options["operating mode"]):
             self.submodels[
                 "leading order external circuit"
             ] = pybamm.external_circuit.LeadingOrderFunctionControl(
-                self.param, self.options["operating mode"]
+                self.param, self.options["operating mode"], self.options
             )
 
     def set_current_collector_submodel(self):
diff --git a/pybamm/models/standard_variables.py b/pybamm/models/standard_variables.py
index b496e1199..0397b34fb 100644
--- a/pybamm/models/standard_variables.py
+++ b/pybamm/models/standard_variables.py
@@ -7,8 +7,9 @@ import numpy as np
 
 class StandardVariables:
     def __init__(self):
-        # Discharge capacity
-        self.Q = pybamm.Variable("Discharge capacity [A.h]")
+        # Discharge capacity and energy
+        self.Q_Ah = pybamm.Variable("Discharge capacity [A.h]")
+        self.Q_Wh = pybamm.Variable("Discharge energy [W.h]")
 
         # Electrolyte concentration
         self.c_e_n = pybamm.Variable(
diff --git a/pybamm/models/submodels/external_circuit/base_external_circuit.py b/pybamm/models/submodels/external_circuit/base_external_circuit.py
index b2116c4cf..e3fe43f0f 100644
--- a/pybamm/models/submodels/external_circuit/base_external_circuit.py
+++ b/pybamm/models/submodels/external_circuit/base_external_circuit.py
@@ -7,32 +7,46 @@ import pybamm
 class BaseModel(pybamm.BaseSubModel):
     """Model to represent the behaviour of the external circuit."""
 
-    def __init__(self, param):
-        super().__init__(param)
+    def __init__(self, param, options):
+        super().__init__(param, options=options)
 
     def get_fundamental_variables(self):
-        Q = pybamm.standard_variables.Q
-        variables = {"Discharge capacity [A.h]": Q}
+        Q_Ah = pybamm.standard_variables.Q_Ah
+        variables = {"Discharge capacity [A.h]": Q_Ah}
+        if self.options["calculate discharge energy"] == "true":
+            Q_Wh = pybamm.standard_variables.Q_Wh
+            variables.update({"Discharge energy [W.h]": Q_Wh})
         return variables
 
     def set_initial_conditions(self, variables):
-        Q = variables["Discharge capacity [A.h]"]
-        self.initial_conditions[Q] = pybamm.Scalar(0)
+        Q_Ah = variables["Discharge capacity [A.h]"]
+        self.initial_conditions[Q_Ah] = pybamm.Scalar(0)
+        if self.options["calculate discharge energy"] == "true":
+            Q_Wh = variables["Discharge energy [W.h]"]
+            self.initial_conditions[Q_Wh] = pybamm.Scalar(0)
 
     def set_rhs(self, variables):
         # ODE for discharge capacity
-        Q = variables["Discharge capacity [A.h]"]
+        Q_Ah = variables["Discharge capacity [A.h]"]
         I = variables["Current [A]"]
-        self.rhs[Q] = I * self.param.timescale / 3600
+
+        self.rhs[Q_Ah] = I * self.param.timescale / 3600
+        if self.options["calculate discharge energy"] == "true":
+            Q_Wh = variables["Discharge energy [W.h]"]
+            V = variables["Terminal voltage [V]"]
+            self.rhs[Q_Wh] = I * V * self.param.timescale / 3600
 
 
 class LeadingOrderBaseModel(BaseModel):
     """Model to represent the behaviour of the external circuit, at leading order."""
 
-    def __init__(self, param):
-        super().__init__(param)
+    def __init__(self, param, options):
+        super().__init__(param, options)
 
     def get_fundamental_variables(self):
-        Q = pybamm.Variable("Leading-order discharge capacity [A.h]")
-        variables = {"Discharge capacity [A.h]": Q}
+        Q_Ah = pybamm.Variable("Leading-order discharge capacity [A.h]")
+        variables = {"Discharge capacity [A.h]": Q_Ah}
+        if self.options["calculate discharge energy"] == "true":
+            Q_Wh = pybamm.Variable("Leading-order discharge energy [W.h]")
+            variables.update({"Discharge energy [W.h]": Q_Wh})
         return variables
diff --git a/pybamm/models/submodels/external_circuit/explicit_control_external_circuit.py b/pybamm/models/submodels/external_circuit/explicit_control_external_circuit.py
index ccf4a0e3f..d2bb8ee63 100644
--- a/pybamm/models/submodels/external_circuit/explicit_control_external_circuit.py
+++ b/pybamm/models/submodels/external_circuit/explicit_control_external_circuit.py
@@ -8,8 +8,8 @@ from .base_external_circuit import BaseModel, LeadingOrderBaseModel
 class ExplicitCurrentControl(BaseModel):
     """External circuit with current control."""
 
-    def __init__(self, param):
-        super().__init__(param)
+    def __init__(self, param, options):
+        super().__init__(param, options)
 
     def get_fundamental_variables(self):
         # Current is given as a function of time
@@ -34,8 +34,8 @@ class ExplicitCurrentControl(BaseModel):
 class ExplicitPowerControl(BaseModel):
     """External circuit with current set explicitly to hit target power."""
 
-    def __init__(self, param):
-        super().__init__(param)
+    def __init__(self, param, options):
+        super().__init__(param, options)
 
     def get_coupled_variables(self, variables):
         param = self.param
@@ -64,8 +64,8 @@ class ExplicitPowerControl(BaseModel):
 class ExplicitResistanceControl(BaseModel):
     """External circuit with current set explicitly to hit target resistance."""
 
-    def __init__(self, param):
-        super().__init__(param)
+    def __init__(self, param, options):
+        super().__init__(param, options)
 
     def get_coupled_variables(self, variables):
         param = self.param
@@ -94,5 +94,5 @@ class ExplicitResistanceControl(BaseModel):
 class LeadingOrderExplicitCurrentControl(ExplicitCurrentControl, LeadingOrderBaseModel):
     """External circuit with current control, for leading order models."""
 
-    def __init__(self, param):
-        super().__init__(param)
+    def __init__(self, param, options):
+        super().__init__(param, options)
diff --git a/pybamm/models/submodels/external_circuit/function_control_external_circuit.py b/pybamm/models/submodels/external_circuit/function_control_external_circuit.py
index 567801603..ad56c9b58 100644
--- a/pybamm/models/submodels/external_circuit/function_control_external_circuit.py
+++ b/pybamm/models/submodels/external_circuit/function_control_external_circuit.py
@@ -16,13 +16,15 @@ class FunctionControl(BaseModel):
         The parameters to use for this submodel
     external_circuit_function : callable
         The function that controls the current
+    options : dict
+        Dictionary of options to use for the submodel
     control : str, optional
         The type of control to use. Must be one of 'algebraic' (default)
         or 'differential'.
     """
 
-    def __init__(self, param, external_circuit_function, control="algebraic"):
-        super().__init__(param)
+    def __init__(self, param, external_circuit_function, options, control="algebraic"):
+        super().__init__(param, options)
         self.external_circuit_function = external_circuit_function
         self.control = control
 
@@ -81,8 +83,8 @@ class VoltageFunctionControl(FunctionControl):
     External circuit with voltage control, implemented as an extra algebraic equation.
     """
 
-    def __init__(self, param):
-        super().__init__(param, self.constant_voltage, control="algebraic")
+    def __init__(self, param, options):
+        super().__init__(param, self.constant_voltage, options, control="algebraic")
 
     def constant_voltage(self, variables):
         V = variables["Terminal voltage [V]"]
@@ -94,8 +96,8 @@ class VoltageFunctionControl(FunctionControl):
 class PowerFunctionControl(FunctionControl):
     """External circuit with power control."""
 
-    def __init__(self, param, control):
-        super().__init__(param, self.constant_power, control=control)
+    def __init__(self, param, options, control):
+        super().__init__(param, self.constant_power, options, control=control)
 
     def constant_power(self, variables):
         I = variables["Current [A]"]
@@ -115,8 +117,8 @@ class PowerFunctionControl(FunctionControl):
 class ResistanceFunctionControl(FunctionControl):
     """External circuit with resistance control."""
 
-    def __init__(self, param, control):
-        super().__init__(param, self.constant_resistance, control=control)
+    def __init__(self, param, options, control):
+        super().__init__(param, self.constant_resistance, options, control=control)
 
     def constant_resistance(self, variables):
         I = variables["Current [A]"]
@@ -146,8 +148,8 @@ class CCCVFunctionControl(FunctionControl):
 
     """
 
-    def __init__(self, param):
-        super().__init__(param, self.cccv, control="differential with max")
+    def __init__(self, param, options):
+        super().__init__(param, self.cccv, options, control="differential with max")
         pybamm.citations.register("Mohtat2021")
 
     def cccv(self, variables):
@@ -165,8 +167,8 @@ class CCCVFunctionControl(FunctionControl):
 class LeadingOrderFunctionControl(FunctionControl, LeadingOrderBaseModel):
     """External circuit with an arbitrary function, at leading order."""
 
-    def __init__(self, param, external_circuit_class, control="algebraic"):
-        super().__init__(param, external_circuit_class, control=control)
+    def __init__(self, param, external_circuit_function, options, control="algebraic"):
+        super().__init__(param, external_circuit_function, options, control=control)
 
     def _get_current_variable(self):
         return pybamm.Variable("Leading-order total current density")
@@ -178,8 +180,8 @@ class LeadingOrderVoltageFunctionControl(LeadingOrderFunctionControl):
     at leading order.
     """
 
-    def __init__(self, param):
-        super().__init__(param, self.constant_voltage, control="algebraic")
+    def __init__(self, param, options):
+        super().__init__(param, self.constant_voltage, options, control="algebraic")
 
     def constant_voltage(self, variables):
         V = variables["Terminal voltage [V]"]
@@ -191,8 +193,8 @@ class LeadingOrderVoltageFunctionControl(LeadingOrderFunctionControl):
 class LeadingOrderPowerFunctionControl(LeadingOrderFunctionControl):
     """External circuit with power control, at leading order."""
 
-    def __init__(self, param):
-        super().__init__(param, self.constant_power, control="algebraic")
+    def __init__(self, param, options):
+        super().__init__(param, self.constant_power, options, control="algebraic")
 
     def constant_power(self, variables):
         I = variables["Current [A]"]
diff --git a/pybamm/simulation.py b/pybamm/simulation.py
index 2f69ec20a..1ea19d103 100644
--- a/pybamm/simulation.py
+++ b/pybamm/simulation.py
@@ -299,7 +299,7 @@ class Simulation:
         # FunctionControl submodel
         # create the FunctionControl submodel and extract variables
         external_circuit_variables = pybamm.external_circuit.FunctionControl(
-            model.param, None
+            model.param, None, model.options
         ).get_fundamental_variables()
 
         # Perform the replacement
@@ -396,7 +396,7 @@ class Simulation:
                     # create the FunctionControl submodel and extract variables
                     external_circuit_variables = (
                         pybamm.external_circuit.FunctionControl(
-                            model.param, None, control=control
+                            model.param, None, model.options, control=control
                         ).get_fundamental_variables()
                     )
 
@@ -470,7 +470,7 @@ class Simulation:
                         new_model.algebraic[
                             i_cell
                         ] = pybamm.external_circuit.VoltageFunctionControl(
-                            new_model.param
+                            new_model.param, model.options
                         ).constant_voltage(
                             new_model.variables
                         )
@@ -478,7 +478,7 @@ class Simulation:
                         new_model.algebraic[
                             i_cell
                         ] = pybamm.external_circuit.PowerFunctionControl(
-                            new_model.param, control="algebraic"
+                            new_model.param, new_model.options, control="algebraic"
                         ).constant_power(
                             new_model.variables
                         )
@@ -486,7 +486,7 @@ class Simulation:
                         new_model.rhs[
                             i_cell
                         ] = pybamm.external_circuit.CCCVFunctionControl(
-                            new_model.param
+                            new_model.param, new_model.options
                         ).cccv(
                             new_model.variables
                         )


Test patch:
diff --git a/tests/unit/test_citations.py b/tests/unit/test_citations.py
index 086147909..34f79ea48 100644
--- a/tests/unit/test_citations.py
+++ b/tests/unit/test_citations.py
@@ -184,7 +184,7 @@ class TestCitations(unittest.TestCase):
 
         citations._reset()
         self.assertNotIn("Mohtat2021", citations._papers_to_cite)
-        pybamm.external_circuit.CCCVFunctionControl(None)
+        pybamm.external_circuit.CCCVFunctionControl(None, None)
         self.assertIn("Mohtat2021", citations._papers_to_cite)
 
     def test_sripad_2020(self):
diff --git a/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py b/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py
index eff8b6490..2a16da398 100644
--- a/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py
+++ b/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py
@@ -14,6 +14,7 @@ OPTIONS_DICT = {
 }
 
 PRINT_OPTIONS_OUTPUT = """\
+'calculate discharge energy': 'false' (possible: ['false', 'true'])
 'cell geometry': 'pouch' (possible: ['arbitrary', 'pouch'])
 'calculate heat source for isothermal models': 'false' (possible: ['false', 'true'])
 'convection': 'none' (possible: ['none', 'uniform transverse', 'full transverse'])
diff --git a/tests/unit/test_models/test_full_battery_models/test_lead_acid/test_loqs.py b/tests/unit/test_models/test_full_battery_models/test_lead_acid/test_loqs.py
index e67d3025e..ba97870ba 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lead_acid/test_loqs.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lead_acid/test_loqs.py
@@ -141,6 +141,11 @@ class TestLeadAcidLOQSExternalCircuits(unittest.TestCase):
         model = pybamm.lead_acid.LOQS(options)
         model.check_well_posedness()
 
+    def test_well_posed_discharge_energy(self):
+        options = {"calculate discharge energy": "true"}
+        model = pybamm.lead_acid.LOQS(options)
+        model.check_well_posedness()
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py
index 1ccc80fd9..b19dcbf80 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py
@@ -219,6 +219,10 @@ class BaseUnitTestLithiumIon:
         }
         self.check_well_posedness(options)
 
+    def test_well_posed_discharge_energy(self):
+        options = {"calculate discharge energy": "true"}
+        self.check_well_posedness(options)
+
     def test_well_posed_external_circuit_voltage(self):
         options = {"operating mode": "voltage"}
         self.check_well_posedness(options)