#871 merge master

CHANGELOG.md

@@ -5,8 +5,12 @@
 -   Added functionality to broadcast to edges ([#891](https://github.com/pybamm-team/PyBaMM/pull/891))
 -   Reformatted and cleaned up `QuickPlot` ([#886](https://github.com/pybamm-team/PyBaMM/pull/886))
 -   Added thermal effects to lead-acid models ([#885](https://github.com/pybamm-team/PyBaMM/pull/885))
+-   Add new symbols `VariableDot`, representing the derivative of a variable wrt time, 
+    and `StateVectorDot`, representing the derivative of a state vector wrt time 
+    ([#858](https://github.com/pybamm-team/PyBaMM/issues/858))
 -   Added additional notebooks showing how to create and compare models ([#877](https://github.com/pybamm-team/PyBaMM/pull/877))
--   Added `Minimum`, `Maximum` and `Sign` operators ([#876](https://github.com/pybamm-team/PyBaMM/pull/876))
+-   Added `Minimum`, `Maximum` and `Sign` operators 
+    ([#876](https://github.com/pybamm-team/PyBaMM/pull/876))
 -   Added a search feature to `FuzzyDict` ([#875](https://github.com/pybamm-team/PyBaMM/pull/875))
 -   Add ambient temperature as a function of time ([#872](https://github.com/pybamm-team/PyBaMM/pull/872))
 -   Added `CasadiAlgebraicSolver` for solving algebraic systems with CasADi ([#868](https://github.com/pybamm-team/PyBaMM/pull/868))

docs/source/expression_tree/state_vector.rst

@@ -3,3 +3,6 @@ State Vector
 
 .. autoclass:: pybamm.StateVector
   :members:
+
+.. autoclass:: pybamm.StateVectorDot
+  :members:

docs/source/expression_tree/variable.rst

@@ -4,5 +4,9 @@ Variable
 .. autoclass:: pybamm.Variable
   :members:
 
+.. autoclass:: pybamm.VariableDot
+  :members:
+
 .. autoclass:: pybamm.ExternalVariable
   :members:
+

examples/notebooks/expression_tree/expression-tree.ipynb

@@ -64,7 +64,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "We can also calculate the expression tree representing the gradient of the equation with respect to $t$ (which is of course simply the scalar value 1),"
+    "We can also calculate the expression tree representing the gradient of the equation with respect to $t$,"
    ]
   },
   {
@@ -84,7 +84,7 @@
     "![](expression_tree2.png)\n",
     "\n",
     "\n",
-    "...and evaluate this expression, which will again give 1."
+    "...and evaluate this expression,"
    ]
   },
   {
@@ -95,7 +95,7 @@
     {
      "data": {
       "text/plain": [
-       "1.0"
+       "array([[-11.]])"
       ]
      },
      "execution_count": 4,
@@ -104,7 +104,7 @@
     }
    ],
    "source": [
-    "diff_wrt_equation.evaluate(1, np.array([2]))"
+    "diff_wrt_equation.evaluate(t=1, y=np.array([2]), y_dot=np.array([2]))"
    ]
   },
   {
@@ -202,6 +202,13 @@
     "\n",
     "After the third stage, our expression tree is now able to be evaluated by one of the solver classes. Note that we have used a single equation above to illustrate the different types of expression trees in PyBaMM, but any given models will consist of many RHS or algebraic equations, along with boundary conditions. See [here](https://github.com/pybamm-team/PyBaMM/blob/master/examples/notebooks/add-model.ipynb) for more details of PyBaMM models."
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {
@@ -220,7 +227,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.3"
+   "version": "3.6.7"
   }
  },
  "nbformat": 4,

examples/scripts/compare_lithium_ion_3D.py

@@ -56,6 +56,6 @@
     solutions[i] = solution
 
 # plot
-output_variables = ["Terminal voltage [V]", "Positive current collector potential"]
+output_variables = ["Terminal voltage [V]", "Negative current collector potential [V]"]
 plot = pybamm.QuickPlot(solutions, output_variables)
 plot.dynamic_plot()

examples/scripts/drive_cycle.py

@@ -3,15 +3,29 @@
 #
 import pybamm
 
+pybamm.set_logging_level("INFO")
+
 # load model and update parameters so the input current is the US06 drive cycle
-model = pybamm.lithium_ion.DFN()
+model = pybamm.lithium_ion.SPMe({"thermal": "x-lumped"})
 param = model.default_parameter_values
 param["Current function [A]"] = "[current data]US06"
 
 # create and run simulation using the CasadiSolver in "fast" mode, remembering to
 # pass in the updated parameters
 sim = pybamm.Simulation(
-    model, parameter_values=param, solver=pybamm.CasadiSolver(mode="fast")
+    model, parameter_values=param, solver=pybamm.CasadiSolver(mode="fast"),
 )
 sim.solve()
-sim.plot()
+sim.plot(
+    [
+        "Negative particle surface concentration [mol.m-3]",
+        "Electrolyte concentration [mol.m-3]",
+        "Positive particle surface concentration [mol.m-3]",
+        "Current [A]",
+        "Negative electrode potential [V]",
+        "Electrolyte potential [V]",
+        "Positive electrode potential [V]",
+        "Terminal voltage [V]",
+        "X-averaged cell temperature",
+    ]
+)

pybamm/__init__.py

@@ -77,11 +77,12 @@ def version(formatted=False):
 from .expression_tree.parameter import Parameter, FunctionParameter
 from .expression_tree.broadcasts import *
 from .expression_tree.scalar import Scalar
-from .expression_tree.variable import Variable, ExternalVariable
+from .expression_tree.variable import Variable, ExternalVariable, VariableDot
+from .expression_tree.variable import VariableBase
 from .expression_tree.independent_variable import *
 from .expression_tree.independent_variable import t
 from .expression_tree.vector import Vector
-from .expression_tree.state_vector import StateVector
+from .expression_tree.state_vector import StateVectorBase, StateVector, StateVectorDot
 
 from .expression_tree.exceptions import *
 

pybamm/discretisations/discretisation.py

@@ -515,6 +515,7 @@ def process_rhs_and_algebraic(self, model):
             equations) and processed_concatenated_algebraic
 
         """
+
         # Discretise right-hand sides, passing domain from variable
         processed_rhs = self.process_dict(model.rhs)
 
@@ -856,6 +857,13 @@ def _process_symbol(self, symbol):
             disc_children = [self.process_symbol(child) for child in symbol.children]
             return symbol._function_new_copy(disc_children)
 
+        elif isinstance(symbol, pybamm.VariableDot):
+            return pybamm.StateVectorDot(
+                *self.y_slices[symbol.get_variable().id],
+                domain=symbol.domain,
+                auxiliary_domains=symbol.auxiliary_domains
+            )
+
         elif isinstance(symbol, pybamm.Variable):
             # Check if variable is a standard variable or an external variable
             if any(symbol.id == var.id for var in self.external_variables.values()):

pybamm/expression_tree/array.py

@@ -98,6 +98,6 @@ def new_copy(self):
             self.entries_string,
         )
 
-    def _base_evaluate(self, t=None, y=None, u=None):
+    def _base_evaluate(self, t=None, y=None, y_dot=None, u=None):
         """ See :meth:`pybamm.Symbol._base_evaluate()`. """
         return self._entries

pybamm/expression_tree/binary_operators.py

@@ -13,7 +13,7 @@ def is_scalar_zero(expr):
     Utility function to test if an expression evaluates to a constant scalar zero
     """
     if expr.is_constant():
-        result = expr.evaluate_ignoring_errors()
+        result = expr.evaluate_ignoring_errors(t=None)
         return isinstance(result, numbers.Number) and result == 0
     else:
         return False
@@ -24,20 +24,20 @@ def is_matrix_zero(expr):
     Utility function to test if an expression evaluates to a constant matrix zero
     """
     if expr.is_constant():
-        result = expr.evaluate_ignoring_errors()
+        result = expr.evaluate_ignoring_errors(t=None)
         return (issparse(result) and result.count_nonzero() == 0) or (
             isinstance(result, np.ndarray) and np.all(result == 0)
         )
     else:
         return False
 
 
-def is_one(expr):
+def is_scalar_one(expr):
     """
     Utility function to test if an expression evaluates to a constant scalar one
     """
     if expr.is_constant():
-        result = expr.evaluate_ignoring_errors()
+        result = expr.evaluate_ignoring_errors(t=None)
         return isinstance(result, numbers.Number) and result == 1
     else:
         return False
@@ -162,21 +162,21 @@ def _binary_new_copy(self, left, right):
         "Default behaviour for new_copy"
         return self.__class__(left, right)
 
-    def evaluate(self, t=None, y=None, u=None, known_evals=None):
+    def evaluate(self, t=None, y=None, y_dot=None, u=None, known_evals=None):
         """ See :meth:`pybamm.Symbol.evaluate()`. """
         if known_evals is not None:
             id = self.id
             try:
                 return known_evals[id], known_evals
             except KeyError:
-                left, known_evals = self.left.evaluate(t, y, u, known_evals)
-                right, known_evals = self.right.evaluate(t, y, u, known_evals)
+                left, known_evals = self.left.evaluate(t, y, y_dot, u, known_evals)
+                right, known_evals = self.right.evaluate(t, y, y_dot, u, known_evals)
                 value = self._binary_evaluate(left, right)
                 known_evals[id] = value
                 return value, known_evals
         else:
-            left = self.left.evaluate(t, y, u)
-            right = self.right.evaluate(t, y, u)
+            left = self.left.evaluate(t, y, y_dot, u)
+            right = self.right.evaluate(t, y, y_dot, u)
             return self._binary_evaluate(left, right)
 
     def _evaluate_for_shape(self):
@@ -252,8 +252,12 @@ def _binary_simplify(self, left, right):
         if is_scalar_zero(right):
             return pybamm.Scalar(1)
 
-        # anything to the power of one is itself
+        # zero to the power of anything is zero
         if is_scalar_zero(left):
+            return pybamm.Scalar(0)
+
+        # anything to the power of one is itself
+        if is_scalar_one(right):
             return left
 
         return self.__class__(left, right)
@@ -425,9 +429,9 @@ def _binary_simplify(self, left, right):
             return zeros_of_shape(shape)
 
         # anything multiplied by a scalar one returns itself
-        if is_one(left):
+        if is_scalar_one(left):
             return right
-        if is_one(right):
+        if is_scalar_one(right):
             return left
 
         return pybamm.simplify_multiplication_division(self.__class__, left, right)
@@ -549,7 +553,7 @@ def _binary_simplify(self, left, right):
                 return pybamm.Array(np.inf * np.ones(left.shape_for_testing))
 
         # anything divided by one is itself
-        if is_one(right):
+        if is_scalar_one(right):
             return left
 
         return pybamm.simplify_multiplication_division(self.__class__, left, right)
@@ -622,9 +626,9 @@ def _binary_simplify(self, left, right):
             return zeros_of_shape(shape)
 
         # anything multiplied by a scalar one returns itself
-        if is_one(left):
+        if is_scalar_one(left):
             return right
-        if is_one(right):
+        if is_scalar_one(right):
             return left
 
         return pybamm.simplify_multiplication_division(self.__class__, left, right)

pybamm/expression_tree/concatenations.py

@@ -54,22 +54,22 @@ def _concatenation_evaluate(self, children_eval):
         else:
             return self.concatenation_function(children_eval)
 
-    def evaluate(self, t=None, y=None, u=None, known_evals=None):
+    def evaluate(self, t=None, y=None, y_dot=None, u=None, known_evals=None):
         """ See :meth:`pybamm.Symbol.evaluate()`. """
         children = self.cached_children
         if known_evals is not None:
             if self.id not in known_evals:
                 children_eval = [None] * len(children)
                 for idx, child in enumerate(children):
                     children_eval[idx], known_evals = child.evaluate(
-                        t, y, u, known_evals
+                        t, y, y_dot, u, known_evals
                     )
                 known_evals[self.id] = self._concatenation_evaluate(children_eval)
             return known_evals[self.id], known_evals
         else:
             children_eval = [None] * len(children)
             for idx, child in enumerate(children):
-                children_eval[idx] = child.evaluate(t, y, u)
+                children_eval[idx] = child.evaluate(t, y, y_dot, u)
             return self._concatenation_evaluate(children_eval)
 
     def new_copy(self):

pybamm/expression_tree/functions.py

@@ -152,19 +152,20 @@ def _function_jac(self, children_jacs):
 
         return jacobian
 
-    def evaluate(self, t=None, y=None, u=None, known_evals=None):
+    def evaluate(self, t=None, y=None, y_dot=None, u=None, known_evals=None):
         """ See :meth:`pybamm.Symbol.evaluate()`. """
         if known_evals is not None:
             if self.id not in known_evals:
                 evaluated_children = [None] * len(self.children)
                 for i, child in enumerate(self.children):
                     evaluated_children[i], known_evals = child.evaluate(
-                        t, y, u, known_evals=known_evals
+                        t, y, y_dot, u, known_evals=known_evals
                     )
                 known_evals[self.id] = self._function_evaluate(evaluated_children)
             return known_evals[self.id], known_evals
         else:
-            evaluated_children = [child.evaluate(t, y, u) for child in self.children]
+            evaluated_children = [child.evaluate(t, y, y_dot, u)
+                                  for child in self.children]
             return self._function_evaluate(evaluated_children)
 
     def _evaluate_for_shape(self):

pybamm/expression_tree/independent_variable.py

@@ -48,7 +48,7 @@ def new_copy(self):
         """ See :meth:`pybamm.Symbol.new_copy()`. """
         return Time()
 
-    def _base_evaluate(self, t, y=None, u=None):
+    def _base_evaluate(self, t=None, y=None, y_dot=None, u=None):
         """ See :meth:`pybamm.Symbol._base_evaluate()`. """
         if t is None:
             raise ValueError("t must be provided")

pybamm/expression_tree/input_parameter.py

@@ -36,7 +36,7 @@ def _jac(self, variable):
         """ See :meth:`pybamm.Symbol._jac()`. """
         return pybamm.Scalar(0)
 
-    def _base_evaluate(self, t=None, y=None, u=None):
+    def _base_evaluate(self, t=None, y=None, y_dot=None, u=None):
         # u should be a dictionary
         # convert 'None' to empty dictionary for more informative error
         if u is None: