Merge pull request #839 from pybamm-team/issue-835-initial-conditions

Issue 835 initial conditions

CHANGELOG.md

@@ -53,7 +53,7 @@
 
 ## Bug fixes
 
--   Time for solver should now be given in seconds ([#832](https://github.com/pybamm-team/PyBaMM/pull/832))
+-   Moved evaluation of initial conditions to solver ([#839](https://github.com/pybamm-team/PyBaMM/pull/839))
 -   Fixed a bug where the first line of the data wasn't loaded when parameters are loaded from data ([#819](https://github.com/pybamm-team/PyBaMM/pull/819))
 -   Made `graphviz` an optional dependency ([#810](https://github.com/pybamm-team/PyBaMM/pull/810))
 -   Fixed examples to run with basic pip installation ([#800](https://github.com/pybamm-team/PyBaMM/pull/800))
@@ -74,6 +74,7 @@
 
 ## Breaking changes
 
+-   Time for solver should now be given in seconds ([#832](https://github.com/pybamm-team/PyBaMM/pull/832))
 -   Model events are now represented as a list of `pybamm.Event` ([#759](https://github.com/pybamm-team/PyBaMM/issues/759)
 -   Removed `ParameterValues.update_model`, whose functionality is now replaced by `InputParameter` ([#801](https://github.com/pybamm-team/PyBaMM/pull/801))
 -   Removed `Outer` and `Kron` nodes as no longer used ([#777](https://github.com/pybamm-team/PyBaMM/pull/777))

pybamm/discretisations/discretisation.py

@@ -178,9 +178,7 @@ def process_model(self, model, inplace=True, check_model=True):
         for event in model.events:
             pybamm.logger.debug("Discretise event '{}'".format(event.name))
             processed_event = pybamm.Event(
-                event.name,
-                self.process_symbol(event.expression),
-                event.event_type
+                event.name, self.process_symbol(event.expression), event.event_type
             )
             processed_events.append(processed_event)
         model_disc.events = processed_events
@@ -402,7 +400,7 @@ def process_initial_conditions(self, model):
         # check that all initial conditions are set
         processed_concatenated_initial_conditions = self._concatenate_in_order(
             processed_initial_conditions, check_complete=True
-        ).evaluate(0, None)
+        )
 
         return processed_initial_conditions, processed_concatenated_initial_conditions
 
@@ -996,17 +994,20 @@ def check_initial_conditions(self, model):
         """Check initial conditions are a numpy array"""
         # Individual
         for var, eqn in model.initial_conditions.items():
-            assert type(eqn.evaluate(0, None)) is np.ndarray, pybamm.ModelError(
+            assert isinstance(
+                eqn.evaluate(t=0, u="shape test"), np.ndarray
+            ), pybamm.ModelError(
                 """
                 initial_conditions must be numpy array after discretisation but they are
                 {} for variable '{}'.
                 """.format(
-                    type(eqn.evaluate(0, None)), var
+                    type(eqn.evaluate(t=0, u="shape test")), var
                 )
             )
         # Concatenated
         assert (
-            type(model.concatenated_initial_conditions) is np.ndarray
+            type(model.concatenated_initial_conditions.evaluate(t=0, u="shape test"))
+            is np.ndarray
         ), pybamm.ModelError(
             """
             Concatenated initial_conditions must be numpy array after discretisation but

pybamm/simulation.py

@@ -474,25 +474,15 @@ def step(
         if solver is None:
             solver = self.solver
 
-        if save is False:
-            # Don't pass previous solution
-            self._solution = solver.step(
-                None,
-                self.built_model,
-                dt,
-                npts=npts,
-                external_variables=external_variables,
-                inputs=inputs,
-            )
-        else:
-            self._solution = solver.step(
-                self._solution,
-                self.built_model,
-                dt,
-                npts=npts,
-                external_variables=external_variables,
-                inputs=inputs,
-            )
+        self._solution = solver.step(
+            self._solution,
+            self.built_model,
+            dt,
+            npts=npts,
+            external_variables=external_variables,
+            inputs=inputs,
+            save=save,
+        )
 
     def get_variable_array(self, *variables):
         """

pybamm/solvers/algebraic_solver.py

@@ -79,7 +79,7 @@ def jacobian(y):
             jacobian = None
 
         # Use "initial conditions" set in model as initial guess
-        y0_guess = model.concatenated_initial_conditions
+        y0_guess = model.concatenated_initial_conditions.evaluate(t=0)
 
         # Solve
         solve_start_time = timer.time()

pybamm/solvers/base_solver.py

@@ -48,7 +48,6 @@ def __init__(
         self.max_steps = max_steps
 
         self.models_set_up = set()
-        self.model_step_times = {}
 
         # Defaults, can be overwritten by specific solver
         self.name = "Base solver"
@@ -115,7 +114,7 @@ def set_up(self, model, inputs=None):
 
         """
         inputs = inputs or {}
-        y0 = model.concatenated_initial_conditions
+        y0 = model.concatenated_initial_conditions.evaluate(0, None, inputs)
 
         # Set model timescale
         model.timescale_eval = model.timescale.evaluate(u=inputs)
@@ -138,9 +137,7 @@ def set_up(self, model, inputs=None):
         if model.convert_to_format != "casadi":
             simp = pybamm.Simplification()
             # Create Jacobian from concatenated rhs and algebraic
-            y = pybamm.StateVector(
-                slice(0, np.size(model.concatenated_initial_conditions))
-            )
+            y = pybamm.StateVector(slice(0, np.size(y0)))
             # set up Jacobian object, for re-use of dict
             jacobian = pybamm.Jacobian()
         else:
@@ -280,13 +277,13 @@ def report(string):
             residuals, residuals_eval, jacobian_eval = process(all_states, "residuals")
             model.residuals_eval = residuals_eval
             model.jacobian_eval = jacobian_eval
-            y0_guess = model.concatenated_initial_conditions.flatten()
+            y0_guess = y0.flatten()
             model.y0 = self.calculate_consistent_state(model, 0, y0_guess)
         else:
             # can use DAE solver to solve ODE model
             model.residuals_eval = Residuals(rhs, "residuals", model)
             model.jacobian_eval = jac_rhs
-            model.y0 = model.concatenated_initial_conditions[:, 0]
+            model.y0 = y0.flatten()
 
         # Save CasADi functions for the CasADi solver
         # Note: when we pass to casadi the ode part of the problem must be in explicit
@@ -564,7 +561,7 @@ def solve(self, model, t_eval, external_variables=None, inputs=None):
 
         # Add model and inputs to solution
         solution.model = model
-        solution.inputs = inputs
+        solution.inputs = ext_and_inputs
 
         # Identify the event that caused termination
         termination = self.get_termination_reason(solution, model.events)
@@ -580,7 +577,14 @@ def solve(self, model, t_eval, external_variables=None, inputs=None):
         return solution
 
     def step(
-        self, old_solution, model, dt, npts=2, external_variables=None, inputs=None
+        self,
+        old_solution,
+        model,
+        dt,
+        npts=2,
+        external_variables=None,
+        inputs=None,
+        save=True,
     ):
         """
         Step the solution of the model forward by a given time increment. The
@@ -604,7 +608,8 @@ def step(
             values at the current time
         inputs : dict, optional
             Any input parameters to pass to the model when solving
-
+        save : bool
+            Turn on to store the solution of all previous timesteps
 
         Raises
         ------
@@ -634,20 +639,22 @@ def step(
         ext_and_inputs = {**external_variables, **inputs}
 
         # Run set up on first step
-        if model not in self.model_step_times:
+        if old_solution is None:
             pybamm.logger.info(
                 "Start stepping {} with {}".format(model.name, self.name)
             )
             self.set_up(model, ext_and_inputs)
-            self.model_step_times[model] = 0.0
+            t = 0.0
             set_up_time = timer.time()
         else:
+            # initialize with old solution
+            t = old_solution.t[-1]
+            model.y0 = old_solution.y[:, -1]
             set_up_time = 0
 
         # Non-dimensionalise dt
         dt_dimensionless = dt / model.timescale_eval
         # Step
-        t = self.model_step_times[model]
         t_eval = np.linspace(t, t + dt_dimensionless, npts)
         # Set inputs and external
         self.set_inputs(model, ext_and_inputs)
@@ -662,15 +669,11 @@ def step(
 
         # Add model and inputs to solution
         solution.model = model
-        solution.inputs = inputs
+        solution.inputs = ext_and_inputs
 
         # Identify the event that caused termination
         termination = self.get_termination_reason(solution, model.events)
 
-        # Set self.t and self.y0 to their values at the final step
-        self.model_step_times[model] = solution.t[-1]
-        model.y0 = solution.y[:, -1]
-
         pybamm.logger.debug("Finish stepping {} ({})".format(model.name, termination))
         if set_up_time:
             pybamm.logger.debug(
@@ -684,7 +687,7 @@ def step(
             pybamm.logger.debug(
                 "Step time: {}".format(timer.format(solution.solve_time))
             )
-        if old_solution is None:
+        if save is False or old_solution is None:
             return solution
         else:
             return old_solution + solution

tests/integration/test_models/standard_model_tests.py

@@ -136,7 +136,7 @@ def evaluate_model(self, simplify=False, use_known_evals=False, to_python=False)
             if simplify:
                 eqn = eqn.simplify()
 
-            y = self.model.concatenated_initial_conditions
+            y = self.model.concatenated_initial_conditions.evaluate(t=0)
             if use_known_evals:
                 eqn_eval, known_evals = eqn.evaluate(0, y, known_evals={})
             elif to_python:

tests/integration/test_solvers/test_external_variables.py

@@ -0,0 +1,65 @@
+#
+# Test solvers with external variables
+#
+import pybamm
+import numpy as np
+import unittest
+
+
+class TestExternalVariables(unittest.TestCase):
+    def test_on_dfn(self):
+        e_height = 0.25
+
+        model = pybamm.lithium_ion.DFN()
+        geometry = model.default_geometry
+        param = model.default_parameter_values
+        param.update({"Electrode height [m]": "[input]"})
+        param.process_model(model)
+        param.process_geometry(geometry)
+        inputs = {"Electrode height [m]": e_height}
+        var = pybamm.standard_spatial_vars
+        var_pts = {var.x_n: 5, var.x_s: 5, var.x_p: 5, var.r_n: 10, var.r_p: 10}
+        spatial_methods = model.default_spatial_methods
+
+        solver = pybamm.CasadiSolver()
+        sim = pybamm.Simulation(
+            model=model,
+            geometry=geometry,
+            parameter_values=param,
+            var_pts=var_pts,
+            spatial_methods=spatial_methods,
+            solver=solver,
+        )
+        sim.solve(t_eval=np.linspace(0, 3600, 100), inputs=inputs)
+
+    def test_external_variables_SPMe(self):
+        model_options = {
+            "thermal": "x-lumped",
+            "external submodels": ["thermal"],
+        }
+        model = pybamm.lithium_ion.SPMe(model_options)
+        sim = pybamm.Simulation(model)
+        t_eval = np.linspace(0, 100, 3)
+        T_av = 0
+        for i in np.arange(1, len(t_eval) - 1):
+            dt = t_eval[i + 1] - t_eval[i]
+            external_variables = {"X-averaged cell temperature": T_av}
+            T_av += 1
+            sim.step(dt, external_variables=external_variables)
+        var = "Terminal voltage [V]"
+        t = sim.solution.t[-1]
+        y = sim.solution.y[:, -1]
+        u = external_variables
+        sim.built_model.variables[var].evaluate(t, y, u)
+        sim.solution[var](t)
+
+
+if __name__ == "__main__":
+    import sys
+
+    print("Add -v for more debug output")
+
+    if "-v" in sys.argv:
+        debug = True
+    pybamm.settings.debug_mode = True
+    unittest.main()

tests/unit/test_discretisations/test_discretisation.py

@@ -593,7 +593,7 @@ def test_process_model_ode(self):
         combined_submesh = mesh.combine_submeshes(*whole_cell)
         disc.process_model(model)
 
-        y0 = model.concatenated_initial_conditions
+        y0 = model.concatenated_initial_conditions.evaluate()
         np.testing.assert_array_equal(
             y0, 3 * np.ones_like(combined_submesh[0].nodes[:, np.newaxis])
         )
@@ -638,7 +638,7 @@ def test_process_model_ode(self):
         model.variables = {"ST": S * T}
 
         disc.process_model(model)
-        y0 = model.concatenated_initial_conditions
+        y0 = model.concatenated_initial_conditions.evaluate()
         y0_expect = np.empty((0, 1))
         for var_id, _ in sorted(disc.y_slices.items(), key=lambda kv: kv[1]):
             if var_id == c.id:
@@ -716,7 +716,7 @@ def test_process_model_dae(self):
         disc.process_model(model)
         combined_submesh = mesh.combine_submeshes(*whole_cell)
 
-        y0 = model.concatenated_initial_conditions
+        y0 = model.concatenated_initial_conditions.evaluate()
         np.testing.assert_array_equal(
             y0,
             np.concatenate(
@@ -815,7 +815,7 @@ def test_process_model_concatenation(self):
         )
 
         disc.process_model(model)
-        y0 = model.concatenated_initial_conditions
+        y0 = model.concatenated_initial_conditions.evaluate()
         np.testing.assert_array_equal(
             y0, 3 * np.ones_like(combined_submesh[0].nodes[:, np.newaxis])
         )
@@ -842,7 +842,7 @@ def test_process_model_not_inplace(self):
         submesh = mesh["negative electrode"]
 
         discretised_model = disc.process_model(model, inplace=False)
-        y0 = discretised_model.concatenated_initial_conditions
+        y0 = discretised_model.concatenated_initial_conditions.evaluate()
         np.testing.assert_array_equal(
             y0, 3 * np.ones_like(submesh[0].nodes[:, np.newaxis])
         )

tests/unit/test_solvers/test_base_solver.py

@@ -41,6 +41,7 @@ def test_ode_solver_fail_with_dae(self):
         model = pybamm.BaseModel()
         a = pybamm.Scalar(1)
         model.algebraic = {a: a}
+        model.concatenated_initial_conditions = pybamm.Scalar(0)
         solver = pybamm.ScipySolver()
         with self.assertRaisesRegex(pybamm.SolverError, "Cannot use ODE solver"):
             solver.set_up(model)

tests/unit/test_solvers/test_scikits_solvers.py

@@ -688,6 +688,7 @@ def test_ode_solver_fail_with_dae(self):
         model = pybamm.BaseModel()
         a = pybamm.Scalar(1)
         model.algebraic = {a: a}
+        model.concatenated_initial_conditions = a
         solver = pybamm.ScikitsOdeSolver()
         with self.assertRaisesRegex(pybamm.SolverError, "Cannot use ODE solver"):
             solver.set_up(model)