Issue 1006 process space var

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## Features
 
+-   Allowed `ProcessedVariable` to handle cases where `len(solution.t)=1` ([#1020](https://github.com/pybamm-team/PyBaMM/pull/1020))
 -   Added `BackwardIndefiniteIntegral` symbol ([#1014](https://github.com/pybamm-team/PyBaMM/pull/1014))
 -   Added `plot` and `plot2D` to enable easy plotting of `pybamm.Array` objects ([#1008](https://github.com/pybamm-team/PyBaMM/pull/1008))
 -   Updated effective current collector models and added example notebook ([#1007](https://github.com/pybamm-team/PyBaMM/pull/1007))
@@ -41,6 +42,7 @@
 
 ## Breaking changes
 
+-   For variables discretised using finite elements the result returned by calling `ProcessedVariable` is now transposed ([#1020](https://github.com/pybamm-team/PyBaMM/pull/1020))
 -   Renamed "surface area density" to "surface area to volume ratio" ([#975](https://github.com/pybamm-team/PyBaMM/pull/975))
 -   Replaced "reaction rate" with "exchange-current density" ([#975](https://github.com/pybamm-team/PyBaMM/pull/975))
 -   Changed the implementation of reactions in submodels ([#948](https://github.com/pybamm-team/PyBaMM/pull/948))

pybamm/solvers/base_solver.py

@@ -616,6 +616,15 @@ def solve(self, model, t_eval=None, external_variables=None, inputs=None):
                 timer.format(solution.total_time),
             )
         )
+
+        # Raise error if solution only contains one timestep (except for algebraic
+        # solvers, where we may only expect one time in the solution)
+        if self.algebraic_solver is False and len(solution.t) == 1:
+            raise pybamm.SolverError(
+                "Solution time vector has length 1. "
+                "Check whether simulation terminated too early."
+            )
+
         return solution
 
     def step(

pybamm/solvers/processed_variable.py

@@ -7,6 +7,26 @@
 import scipy.interpolate as interp
 
 
+def make_interp2D_fun(input, interpolant):
+    """
+    Calls and returns a 2D interpolant of the correct shape depending on the
+    shape of the input
+    """
+    first_dim, second_dim, _ = input
+    if isinstance(first_dim, np.ndarray) and isinstance(second_dim, np.ndarray):
+        first_dim = first_dim[:, 0, 0]
+        second_dim = second_dim[:, 0]
+        return interpolant(second_dim, first_dim)
+    elif isinstance(first_dim, np.ndarray):
+        first_dim = first_dim[:, 0]
+        return interpolant(second_dim, first_dim)[:, 0]
+    elif isinstance(second_dim, np.ndarray):
+        second_dim = second_dim[:, 0]
+        return interpolant(second_dim, first_dim)
+    else:
+        return interpolant(second_dim, first_dim)[0]
+
+
 class ProcessedVariable(object):
     """
     An object that can be evaluated at arbitrary (scalars or vectors) t and x, and
@@ -55,48 +75,11 @@ def __init__(self, base_variable, solution, known_evals=None):
             and "current collector" in self.domain
             and isinstance(self.mesh[0], pybamm.ScikitSubMesh2D)
         ):
-            if len(solution.t) == 1:
-                # space only (steady solution)
-                self.initialise_2D_fixed_t_scikit_fem()
-            else:
-                self.initialise_2D_scikit_fem()
+            self.initialise_2D_scikit_fem()
 
         # check variable shape
         else:
-            if len(solution.t) == 1:
-                # Implementing a workaround for 0D and 1D variables. Processing
-                # variables that are functions of space only needs to be implemented
-                # properly, see #1006
-                if (
-                    isinstance(self.base_eval, numbers.Number)
-                    or len(self.base_eval.shape) == 0
-                    or self.base_eval.shape[0] == 1
-                ):
-                    # Scalar value
-                    t = self.t_sol
-                    u = self.u_sol
-                    inputs = {name: inp[0] for name, inp in self.inputs.items()}
-
-                    entries = self.base_variable.evaluate(t, u, inputs=inputs)
-
-                    def fun(t):
-                        return entries
-
-                    self._interpolation_function = fun
-                    self.entries = entries
-                    self.dimensions = 0
-                else:
-                    # 1D function of space only
-                    n = self.mesh[0].npts
-                    base_shape = self.base_eval.shape[0]
-                    if base_shape in [n, n + 1]:
-                        self.initialise_1D(fixed_t=True)
-                    else:
-                        raise pybamm.SolverError(
-                            "Solution time vector must have length > 1. Check whether "
-                            "simulation terminated too early."
-                        )
-            elif (
+            if (
                 isinstance(self.base_eval, numbers.Number)
                 or len(self.base_eval.shape) == 0
                 or self.base_eval.shape[0] == 1
@@ -141,9 +124,21 @@ def initialise_0D(self):
                 entries[idx] = self.base_variable.evaluate(t, u, inputs=inputs)
 
         # No discretisation provided, or variable has no domain (function of t only)
-        self._interpolation_function = interp.interp1d(
-            self.t_sol, entries, kind="linear", fill_value=np.nan, bounds_error=False
-        )
+        if len(self.t_sol) == 1:
+            # Variable is just a scalar value, but we need to create a callable
+            # function to be consitent with other processed variables
+            def fun(t):
+                return entries
+
+            self._interpolation_function = fun
+        else:
+            self._interpolation_function = interp.interp1d(
+                self.t_sol,
+                entries,
+                kind="linear",
+                fill_value=np.nan,
+                bounds_error=False,
+            )
 
         self.entries = entries
         self.dimensions = 0
@@ -208,19 +203,22 @@ def initialise_1D(self, fixed_t=False):
         self.internal_boundaries = self.mesh[0].internal_boundaries
 
         # set up interpolation
-        # note that the order of 't' and 'space' is the reverse of what you'd expect
-        # TODO: fix processing when variable is a function of space only
-        if fixed_t:
-            # Hack for 1D space only
+        if len(self.t_sol) == 1:
+            # function of space only
             interpolant = interp.interp1d(
                 space, entries_for_interp[:, 0], kind="linear", fill_value=np.nan
             )
 
             def interp_fun(t, z):
-                return interpolant(z)[:, np.newaxis]
+                if isinstance(z, np.ndarray):
+                    return interpolant(z)[:, np.newaxis]
+                else:
+                    return interpolant(z)
 
             self._interpolation_function = interp_fun
         else:
+            # function of space and time. Note that the order of 't' and 'space'
+            # is the reverse of what you'd expect
             self._interpolation_function = interp.interp2d(
                 self.t_sol, space, entries_for_interp, kind="linear", fill_value=np.nan
             )
@@ -297,40 +295,29 @@ def initialise_2D(self):
         self.second_dim_pts = second_dim_pts
 
         # set up interpolation
-        self._interpolation_function = interp.RegularGridInterpolator(
-            (first_dim_pts, second_dim_pts, self.t_sol),
-            entries,
-            method="linear",
-            fill_value=np.nan,
-        )
-
-    def initialise_2D_fixed_t_scikit_fem(self):
-        y_sol = self.mesh[0].edges["y"]
-        len_y = len(y_sol)
-        z_sol = self.mesh[0].edges["z"]
-        len_z = len(z_sol)
-
-        # Evaluate the base_variable
-        inputs = {name: inp[0] for name, inp in self.inputs.items()}
+        if len(self.t_sol) == 1:
+            # function of space only. Note the order of the points is the reverse
+            # of what you'd expect
+            interpolant = interp.interp2d(
+                second_dim_pts,
+                first_dim_pts,
+                entries[:, :, 0],
+                kind="linear",
+                fill_value=np.nan,
+            )
 
-        entries = np.reshape(
-            self.base_variable.evaluate(0, self.u_sol, inputs=inputs), [len_y, len_z]
-        )
+            def interp_fun(input):
+                return make_interp2D_fun(input, interpolant)
 
-        # assign attributes for reference
-        self.entries = entries
-        self.dimensions = 2
-        self.y_sol = y_sol
-        self.z_sol = z_sol
-        self.first_dimension = "y"
-        self.second_dimension = "z"
-        self.first_dim_pts = y_sol
-        self.second_dim_pts = z_sol
-
-        # set up interpolation
-        self._interpolation_function = interp.interp2d(
-            y_sol, z_sol, entries, kind="linear", fill_value=np.nan
-        )
+            self._interpolation_function = interp_fun
+        else:
+            # function of space and time.
+            self._interpolation_function = interp.RegularGridInterpolator(
+                (first_dim_pts, second_dim_pts, self.t_sol),
+                entries,
+                method="linear",
+                fill_value=np.nan,
+            )
 
     def initialise_2D_scikit_fem(self):
         y_sol = self.mesh[0].edges["y"]
@@ -349,11 +336,15 @@ def initialise_2D_scikit_fem(self):
                 eval_and_known_evals = self.base_variable.evaluate(
                     t, u, inputs=inputs, known_evals=self.known_evals[t]
                 )
-                entries[:, :, idx] = np.reshape(eval_and_known_evals[0], [len_y, len_z])
+                entries[:, :, idx] = np.reshape(
+                    eval_and_known_evals[0], [len_y, len_z], order="F"
+                )
                 self.known_evals[t] = eval_and_known_evals[1]
             else:
                 entries[:, :, idx] = np.reshape(
-                    self.base_variable.evaluate(t, u, inputs=inputs), [len_y, len_z]
+                    self.base_variable.evaluate(t, u, inputs=inputs),
+                    [len_y, len_z],
+                    order="F",
                 )
 
         # assign attributes for reference
@@ -367,23 +358,49 @@ def initialise_2D_scikit_fem(self):
         self.second_dim_pts = z_sol
 
         # set up interpolation
-        self._interpolation_function = interp.RegularGridInterpolator(
-            (y_sol, z_sol, self.t_sol), entries, method="linear", fill_value=np.nan
-        )
+        if len(self.t_sol) == 1:
+            # function of space only. Note the order of the points is the reverse
+            # of what you'd expect
+            interpolant = interp.interp2d(
+                z_sol, y_sol, entries, kind="linear", fill_value=np.nan
+            )
+
+            def interp_fun(input):
+                return make_interp2D_fun(input, interpolant)
+
+            self._interpolation_function = interp_fun
+        else:
+            # function of space and time.
+            self._interpolation_function = interp.RegularGridInterpolator(
+                (y_sol, z_sol, self.t_sol), entries, method="linear", fill_value=np.nan
+            )
 
     def __call__(self, t=None, x=None, r=None, y=None, z=None, warn=True):
         """
         Evaluate the variable at arbitrary t (and x, r, y and/or z), using interpolation
         """
+        # If t is None and there is only one value of time in the soluton (i.e.
+        # the solution is independent of time) then we set t equal to the value
+        # stored in the solution. If the variable is constant (doesn't depend on
+        # time) evaluate arbitrarily at the first value of t. Otherwise, raise
+        # an error
+        if t is None:
+            if len(self.t_sol) == 1:
+                t = self.t_sol
+            elif self.base_variable.is_constant():
+                t = self.t_sol[0]
+            else:
+                raise ValueError(
+                    "t cannot be None for variable {}".format(self.base_variable)
+                )
+
+        # Call interpolant of correct spatial dimension
         if self.dimensions == 0:
             out = self._interpolation_function(t)
         elif self.dimensions == 1:
             out = self.call_1D(t, x, r, z)
         elif self.dimensions == 2:
-            if t is None:
-                out = self._interpolation_function(y, z)
-            else:
-                out = self.call_2D(t, x, r, y, z)
+            out = self.call_2D(t, x, r, y, z)
         if warn is True and np.isnan(out).any():
             pybamm.logger.warning(
                 "Calling variable outside interpolation range (returns 'nan')"

tests/unit/test_solvers/test_base_solver.py

@@ -71,6 +71,17 @@ def test_nonmonotonic_teval(self):
         ):
             solver.solve(model, np.array([1, 2, 3, 2]))
 
+    def test_solution_time_length_fail(self):
+        model = pybamm.BaseModel()
+        v = pybamm.Scalar(1)
+        model.variables = {"v": v}
+        solver = pybamm.DummySolver()
+        t_eval = np.array([0])
+        with self.assertRaisesRegex(
+            pybamm.SolverError, "Solution time vector has length 1"
+        ):
+            solver.solve(model, t_eval)
+
     def test_block_symbolic_inputs(self):
         solver = pybamm.BaseSolver(rtol=1e-2, atol=1e-4)
         model = pybamm.BaseModel()
@@ -203,13 +214,13 @@ def algebraic_eval(self, t, y, inputs):
             solver.calculate_consistent_state(Model())
         solver = pybamm.BaseSolver(root_method="lm")
         with self.assertRaisesRegex(
-            pybamm.SolverError, "Could not find acceptable solution: solver terminated",
+            pybamm.SolverError, "Could not find acceptable solution: solver terminated"
         ):
             solver.calculate_consistent_state(Model())
         # with casadi
         solver = pybamm.BaseSolver(root_method="casadi")
         with self.assertRaisesRegex(
-            pybamm.SolverError, "Could not find acceptable solution: .../casadi",
+            pybamm.SolverError, "Could not find acceptable solution: .../casadi"
         ):
             solver.calculate_consistent_state(Model())