Issue 1015 algebraic solver

CHANGELOG.md

@@ -4,6 +4,7 @@
 
 -   Reformatted Landesfeind electrolytes ([#1064](https://github.com/pybamm-team/PyBaMM/pull/1064))
 -   Adapted examples to be run in Google Colab ([#1061](https://github.com/pybamm-team/PyBaMM/pull/1061))
+-   Added some new solvers for algebraic models ([#1059](https://github.com/pybamm-team/PyBaMM/pull/1059))
 -   Added `length_scales` attribute to models ([#1058](https://github.com/pybamm-team/PyBaMM/pull/1058))
 -   Added averaging in secondary dimensions ([#1057](https://github.com/pybamm-team/PyBaMM/pull/1057))
 

pybamm/discretisations/discretisation.py

@@ -160,8 +160,6 @@ def process_model(self, model, inplace=True, check_model=True):
         # Set the y split for variables
         pybamm.logger.info("Set variable slices for {}".format(model.name))
         self.set_variable_slices(variables)
-        # Keep a record of y_slices in the model
-        model.y_slices = self.y_slices_explicit
 
         # now add extrapolated external variables to the boundary conditions
         # if required by the spatial method
@@ -183,6 +181,11 @@ def process_model(self, model, inplace=True, check_model=True):
             # create an empty copy of the original model
             model_disc = model.new_copy()
 
+        # Keep a record of y_slices in the model
+        model_disc.y_slices = self.y_slices_explicit
+        # Keep a record of the bounds in the model
+        model_disc.bounds = self.bounds
+
         model_disc.bcs = self.bcs
 
         pybamm.logger.info("Discretise initial conditions for {}".format(model.name))
@@ -240,11 +243,13 @@ def set_variable_slices(self, variables):
         variables : iterable of :class:`pybamm.Variables`
             The variables for which to set slices
         """
-        # Set up y_slices
+        # Set up y_slices and bounds
         y_slices = defaultdict(list)
         y_slices_explicit = defaultdict(list)
         start = 0
         end = 0
+        lower_bounds = []
+        upper_bounds = []
         # Iterate through unpacked variables, adding appropriate slices to y_slices
         for variable in variables:
             # Add up the size of all the domains in variable.domain
@@ -261,20 +266,33 @@ def set_variable_slices(self, variables):
                     for child, mesh in meshes.items():
                         for domain_mesh in mesh:
                             end += domain_mesh.npts_for_broadcast_to_nodes
+                        # Add to slices
                         y_slices[child.id].append(slice(start, end))
                         y_slices_explicit[child].append(slice(start, end))
+                        # Add to bounds
+                        lower_bounds.extend([child.bounds[0]] * (end - start))
+                        upper_bounds.extend([child.bounds[1]] * (end - start))
+                        # Increment start
                         start = end
             else:
                 end += self._get_variable_size(variable)
+                # Add to slices
                 y_slices[variable.id].append(slice(start, end))
                 y_slices_explicit[variable].append(slice(start, end))
+                # Add to bounds
+                lower_bounds.extend([variable.bounds[0]] * (end - start))
+                upper_bounds.extend([variable.bounds[1]] * (end - start))
+                # Increment start
                 start = end
 
         # Convert y_slices back to normal dictionary
         self.y_slices = dict(y_slices)
         # Also keep a record of what the y_slices are, to be stored in the model
         self.y_slices_explicit = dict(y_slices_explicit)
 
+        # Also keep a record of bounds
+        self.bounds = (np.array(lower_bounds), np.array(upper_bounds))
+
         # reset discretised_symbols
         self._discretised_symbols = {}
 
@@ -885,7 +903,7 @@ def _process_symbol(self, symbol):
                 # Broadcast new_child to the domain specified by symbol.domain
                 # Different discretisations may broadcast differently
                 if symbol.domain == []:
-                    symbol = disc_child * pybamm.Vector(np.array([1]))
+                    symbol = disc_child * pybamm.Vector([1])
                 else:
                     symbol = spatial_method.broadcast(
                         disc_child,
@@ -921,7 +939,7 @@ def _process_symbol(self, symbol):
             return pybamm.StateVectorDot(
                 *self.y_slices[symbol.get_variable().id],
                 domain=symbol.domain,
-                auxiliary_domains=symbol.auxiliary_domains
+                auxiliary_domains=symbol.auxiliary_domains,
             )
 
         elif isinstance(symbol, pybamm.Variable):
@@ -972,7 +990,7 @@ def _process_symbol(self, symbol):
                 return pybamm.StateVector(
                     *y_slices,
                     domain=symbol.domain,
-                    auxiliary_domains=symbol.auxiliary_domains
+                    auxiliary_domains=symbol.auxiliary_domains,
                 )
 
         elif isinstance(symbol, pybamm.SpatialVariable):

pybamm/expression_tree/array.py

@@ -13,8 +13,9 @@ class Array(pybamm.Symbol):
     Parameters
     ----------
 
-    entries : numpy.array
-        the array associated with the node
+    entries : numpy.array or list
+        the array associated with the node. If a list is provided, it is converted to a
+        numpy array
     name : str, optional
         the name of the node
     domain : iterable of str, optional
@@ -35,6 +36,8 @@ def __init__(
         auxiliary_domains=None,
         entries_string=None,
     ):
+        if isinstance(entries, list):
+            entries = np.array(entries)
         if entries.ndim == 1:
             entries = entries[:, np.newaxis]
         if name is None:

pybamm/expression_tree/concatenations.py

@@ -127,7 +127,7 @@ def __init__(self, *children):
         # so that we can concatenate them
         for i, child in enumerate(children):
             if child.evaluates_to_number():
-                children[i] = child * pybamm.Vector(np.array([1]))
+                children[i] = child * pybamm.Vector([1])
         super().__init__(
             *children,
             name="numpy concatenation",

pybamm/expression_tree/matrix.py

@@ -2,7 +2,7 @@
 # Matrix class
 #
 import pybamm
-
+import numpy as np
 from scipy.sparse import issparse
 
 
@@ -21,6 +21,8 @@ def __init__(
         auxiliary_domains=None,
         entries_string=None,
     ):
+        if isinstance(entries, list):
+            entries = np.array(entries)
         if name is None:
             name = "Matrix {!s}".format(entries.shape)
             if issparse(entries):

pybamm/expression_tree/state_vector.py

@@ -26,7 +26,7 @@ class StateVectorBase(pybamm.Symbol):
         List of boolean arrays representing slices. Default is None, in which case the
         evaluation_array is computed from y_slices.
 
-    *Extends:* :class:`Array`
+    *Extends:* :class:`pybamm.Symbol`
     """
 
     def __init__(
@@ -212,7 +212,7 @@ class StateVector(StateVectorBase):
         List of boolean arrays representing slices. Default is None, in which case the
         evaluation_array is computed from y_slices.
 
-    *Extends:* :class:`Array`
+    *Extends:* :class:`pybamm.StateVectorBase`
     """
 
     def __init__(

pybamm/expression_tree/variable.py

@@ -26,20 +26,33 @@ class VariableBase(pybamm.Symbol):
         collector'. For the DFN, the particle concentration would be a Variable with
         domain 'negative particle', secondary domain 'negative electrode' and tertiary
         domain 'current collector'
+    bounds : tuple, optional
+        Physical bounds on the variable
 
     *Extends:* :class:`Symbol`
     """
 
-    def __init__(self, name, domain=None, auxiliary_domains=None):
+    def __init__(self, name, domain=None, auxiliary_domains=None, bounds=None):
         if domain is None:
             domain = []
         if auxiliary_domains is None:
             auxiliary_domains = {}
         super().__init__(name, domain=domain, auxiliary_domains=auxiliary_domains)
+        if bounds is None:
+            bounds = (-np.inf, np.inf)
+        else:
+            if bounds[0] >= bounds[1]:
+                raise ValueError(
+                    "Invalid bounds {}. ".format(bounds)
+                    + "Lower bound should be strictly less than upper bound."
+                )
+        self.bounds = bounds
 
     def new_copy(self):
         """ See :meth:`pybamm.Symbol.new_copy()`. """
-        return self.__class__(self.name, self.domain, self.auxiliary_domains)
+        return self.__class__(
+            self.name, self.domain, self.auxiliary_domains, self.bounds
+        )
 
     def _evaluate_for_shape(self):
         """ See :meth:`pybamm.Symbol.evaluate_for_shape_using_domain()` """
@@ -59,21 +72,24 @@ class Variable(VariableBase):
 
     name : str
         name of the node
-    domain : iterable of str
+        domain : iterable of str, optional
         list of domains that this variable is valid over
-    auxiliary_domains : dict
+    auxiliary_domains : dict, optional
         dictionary of auxiliary domains ({'secondary': ..., 'tertiary': ...}). For
         example, for the single particle model, the particle concentration would be a
         Variable with domain 'negative particle' and secondary auxiliary domain 'current
         collector'. For the DFN, the particle concentration would be a Variable with
         domain 'negative particle', secondary domain 'negative electrode' and tertiary
         domain 'current collector'
-
+    bounds : tuple, optional
+        Physical bounds on the variable
     *Extends:* :class:`Symbol`
     """
 
-    def __init__(self, name, domain=None, auxiliary_domains=None):
-        super().__init__(name, domain=domain, auxiliary_domains=auxiliary_domains)
+    def __init__(self, name, domain=None, auxiliary_domains=None, bounds=None):
+        super().__init__(
+            name, domain=domain, auxiliary_domains=auxiliary_domains, bounds=bounds
+        )
 
     def diff(self, variable):
         if variable.id == self.id:
@@ -110,11 +126,13 @@ class VariableDot(VariableBase):
         collector'. For the DFN, the particle concentration would be a Variable with
         domain 'negative particle', secondary domain 'negative electrode' and tertiary
         domain 'current collector'
-
+    bounds : tuple, optional
+        Physical bounds on the variable. Included for compatibility with `VariableBase`,
+        but ignored.
     *Extends:* :class:`Symbol`
     """
 
-    def __init__(self, name, domain=None, auxiliary_domains=None):
+    def __init__(self, name, domain=None, auxiliary_domains=None, bounds=None):
         super().__init__(name, domain=domain, auxiliary_domains=auxiliary_domains)
 
     def get_variable(self):

pybamm/expression_tree/vector.py

@@ -21,6 +21,8 @@ def __init__(
         auxiliary_domains=None,
         entries_string=None,
     ):
+        if isinstance(entries, list):
+            entries = np.array(entries)
         # make sure that entries are a vector (can be a column vector)
         if entries.ndim == 1:
             entries = entries[:, np.newaxis]

pybamm/models/standard_variables.py

@@ -2,27 +2,33 @@
 # Standard variables for the models
 #
 import pybamm
+import numpy as np
 
 # Electrolyte concentration
 c_e_n = pybamm.Variable(
     "Negative electrolyte concentration",
     domain="negative electrode",
     auxiliary_domains={"secondary": "current collector"},
+    bounds=(0, np.inf),
 )
 c_e_s = pybamm.Variable(
     "Separator electrolyte concentration",
     domain="separator",
     auxiliary_domains={"secondary": "current collector"},
+    bounds=(0, np.inf),
 )
 c_e_p = pybamm.Variable(
     "Positive electrolyte concentration",
     domain="positive electrode",
     auxiliary_domains={"secondary": "current collector"},
+    bounds=(0, np.inf),
 )
 c_e = pybamm.Concatenation(c_e_n, c_e_s, c_e_p)
 
 c_e_av = pybamm.Variable(
-    "X-averaged electrolyte concentration", domain="current collector"
+    "X-averaged electrolyte concentration",
+    domain="current collector",
+    bounds=(0, np.inf),
 )
 
 # Electrolyte potential
@@ -105,6 +111,7 @@
         "secondary": "negative electrode",
         "tertiary": "current collector",
     },
+    bounds=(0, 1),
 )
 c_s_p = pybamm.Variable(
     "Positive particle concentration",
@@ -113,32 +120,41 @@
         "secondary": "positive electrode",
         "tertiary": "current collector",
     },
+    bounds=(0, 1),
 )
 c_s_n_xav = pybamm.Variable(
     "X-averaged negative particle concentration",
     domain="negative particle",
     auxiliary_domains={"secondary": "current collector"},
+    bounds=(0, 1),
 )
 c_s_p_xav = pybamm.Variable(
     "X-averaged positive particle concentration",
     domain="positive particle",
     auxiliary_domains={"secondary": "current collector"},
+    bounds=(0, 1),
 )
 c_s_n_surf = pybamm.Variable(
     "Negative particle surface concentration",
     domain="negative electrode",
     auxiliary_domains={"secondary": "current collector"},
+    bounds=(0, 1),
 )
 c_s_p_surf = pybamm.Variable(
     "Positive particle surface concentration",
     domain="positive electrode",
     auxiliary_domains={"secondary": "current collector"},
+    bounds=(0, 1),
 )
 c_s_n_surf_xav = pybamm.Variable(
-    "X-averaged negative particle surface concentration", domain="current collector"
+    "X-averaged negative particle surface concentration",
+    domain="current collector",
+    bounds=(0, 1),
 )
 c_s_p_surf_xav = pybamm.Variable(
-    "X-averaged positive particle surface concentration", domain="current collector"
+    "X-averaged positive particle surface concentration",
+    domain="current collector",
+    bounds=(0, 1),
 )
 
 
@@ -147,26 +163,31 @@
     "Negative electrode porosity",
     domain="negative electrode",
     auxiliary_domains={"secondary": "current collector"},
+    bounds=(0, 1),
 )
 eps_s = pybamm.Variable(
     "Separator porosity",
     domain="separator",
     auxiliary_domains={"secondary": "current collector"},
+    bounds=(0, 1),
 )
 eps_p = pybamm.Variable(
     "Positive electrode porosity",
     domain="positive electrode",
     auxiliary_domains={"secondary": "current collector"},
+    bounds=(0, 1),
 )
 eps = pybamm.Concatenation(eps_n, eps_s, eps_p)
 
 # Piecewise constant (for asymptotic models)
 eps_n_pc = pybamm.Variable(
-    "X-averaged negative electrode porosity", domain="current collector"
+    "X-averaged negative electrode porosity", domain="current collector", bounds=(0, 1),
+)
+eps_s_pc = pybamm.Variable(
+    "X-averaged separator porosity", domain="current collector", bounds=(0, 1)
 )
-eps_s_pc = pybamm.Variable("X-averaged separator porosity", domain="current collector")
 eps_p_pc = pybamm.Variable(
-    "X-averaged positive electrode porosity", domain="current collector"
+    "X-averaged positive electrode porosity", domain="current collector", bounds=(0, 1),
 )
 
 eps_piecewise_constant = pybamm.Concatenation(
@@ -219,4 +240,3 @@
     domain=["negative electrode"],
     auxiliary_domains={"secondary": "current collector"},
 )
-