#775 get external test working

Author: valentinsulzer

pybamm/solvers/base_solver.py

@@ -49,7 +49,7 @@ def atol(self):
     def atol(self, value):
         self._atol = value
 
-    def solve(self, model, t_eval, inputs=None):
+    def solve(self, model, t_eval, external_variables=None, inputs=None):
         """
         Execute the solver setup and calculate the solution of the model at
         specified times.
@@ -61,6 +61,9 @@ def solve(self, model, t_eval, inputs=None):
             initial_conditions
         t_eval : numeric type
             The times at which to compute the solution
+        external_variables : dict
+            A dictionary of external variables and their corresponding
+            values at the current time
         inputs : dict, optional
             Any input parameters to pass to the model when solving
 
@@ -80,6 +83,8 @@ def solve(self, model, t_eval, inputs=None):
         timer = pybamm.Timer()
         start_time = timer.time()
         inputs = inputs or {}
+        self.y_pad = np.zeros((model.y_length - model.external_start, 1))
+        self.set_external_variables(model, external_variables)
         if model.convert_to_format == "casadi" or isinstance(self, pybamm.CasadiSolver):
             self.set_up_casadi(model, inputs)
         else:

pybamm/solvers/casadi_solver.py

@@ -68,7 +68,7 @@ def __init__(
         self.extra_options = extra_options
         self.name = "CasADi solver ({}) with '{}' mode".format(method, mode)
 
-    def solve(self, model, t_eval, inputs=None):
+    def solve(self, model, t_eval, external_variables=None, inputs=None):
         """
         Execute the solver setup and calculate the solution of the model at
         specified times.
@@ -80,6 +80,9 @@ def solve(self, model, t_eval, inputs=None):
             initial_conditions
         t_eval : numeric type
             The times at which to compute the solution
+        external_variables : dict
+            A dictionary of external variables and their corresponding
+            values at the current time
         inputs : dict, optional
             Any input parameters to pass to the model when solving
 
@@ -93,11 +96,15 @@ def solve(self, model, t_eval, inputs=None):
         """
         if self.mode == "fast":
             # Solve model normally by calling the solve method from parent class
-            return super().solve(model, t_eval, inputs=inputs)
+            return super().solve(
+                model, t_eval, external_variables=external_variables, inputs=inputs
+            )
         elif model.events == {}:
             pybamm.logger.info("No events found, running fast mode")
             # Solve model normally by calling the solve method from parent class
-            return super().solve(model, t_eval, inputs=inputs)
+            return super().solve(
+                model, t_eval, external_variables=external_variables, inputs=inputs
+            )
         elif self.mode == "safe":
             # Step-and-check
             timer = pybamm.Timer()
@@ -122,7 +129,12 @@ def solve(self, model, t_eval, inputs=None):
                     # different to t_eval, but shouldn't matter too much as it should
                     # only happen near events.
                     try:
-                        current_step_sol = self.step(model, dt, inputs=inputs)
+                        current_step_sol = self.step(
+                            model,
+                            dt,
+                            external_variables=external_variables,
+                            inputs=inputs,
+                        )
                         solved = True
                     except pybamm.SolverError:
                         dt /= 2
@@ -229,6 +241,11 @@ def integrate_casadi(self, rhs, algebraic, y0, t_eval, inputs=None):
             Any input parameters to pass to the model when solving
         """
         inputs = inputs or {}
+        if self.y_ext is None:
+            y_ext = np.array([])
+        else:
+            y_ext = self.y_ext
+
         options = {
             "grid": t_eval,
             "reltol": self.rtol,
@@ -242,16 +259,15 @@ def integrate_casadi(self, rhs, algebraic, y0, t_eval, inputs=None):
         # set up and solve
         t = casadi.MX.sym("t")
         u = casadi.vertcat(*[x for x in inputs.values()])
-        y_diff = self.y_diff
+        y0_w_ext = casadi.vertcat(y0, y_ext[len(y0) :])
+        y_diff = casadi.MX.sym("y_diff", rhs(0, y0_w_ext, u).shape[0])
         problem = {"t": t, "x": y_diff}
         if algebraic is None:
-            y_casadi_w_ext = casadi.vertcat(y_diff, self.y_ext[y_diff.shape[0] :])
+            y_casadi_w_ext = casadi.vertcat(y_diff, y_ext[len(y0) :])
             problem.update({"ode": rhs(t, y_casadi_w_ext, u)})
         else:
             y_alg = self.y_alg
-            y_casadi_w_ext = casadi.vertcat(
-                y_diff, y_alg, self.y_ext[y_diff.shape[0] + y_alg.shape[0] :]
-            )
+            y_casadi_w_ext = casadi.vertcat(y_diff, y_alg, y_ext[len(y0) :])
             problem.update(
                 {
                     "z": y_alg,

tests/unit/test_solvers/test_casadi_solver.py

@@ -214,6 +214,30 @@ def test_model_solver_with_inputs(self):
         np.testing.assert_array_equal(solution.t, t_eval[: len(solution.t)])
         np.testing.assert_allclose(solution.y[0], np.exp(-0.1 * solution.t), rtol=1e-06)
 
+    def test_model_solver_with_external(self):
+        # Create model
+        model = pybamm.BaseModel()
+        domain = ["negative electrode", "separator", "positive electrode"]
+        var1 = pybamm.Variable("var1", domain=domain)
+        var2 = pybamm.Variable("var2", domain=domain)
+        model.rhs = {var1: -var2}
+        model.initial_conditions = {var1: 1}
+        model.external_variables = [var2]
+        model.variables = {"var1": var1, "var2": var2}
+        # No need to set parameters; can use base discretisation (no spatial
+        # operators)
+
+        # create discretisation
+        mesh = get_mesh_for_testing()
+        spatial_methods = {"macroscale": pybamm.FiniteVolume()}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+        disc.process_model(model)
+        # Solve
+        solver = pybamm.CasadiSolver(rtol=1e-8, atol=1e-8)
+        t_eval = np.linspace(0, 10, 100)
+        solution = solver.solve(model, t_eval, external_variables={"var2": 0.5})
+        np.testing.assert_allclose(solution.y[0], 1 - 0.5 * solution.t, rtol=1e-06)
+
     def test_model_solver_with_non_identity_mass(self):
         model = pybamm.BaseModel()
         var1 = pybamm.Variable("var1", domain="negative electrode")

tests/unit/test_solvers/test_scikits_solvers.py

@@ -745,6 +745,30 @@ def test_model_solver_dae_inputs_events(self):
             np.testing.assert_allclose(solution.y[0], np.exp(0.1 * solution.t))
             np.testing.assert_allclose(solution.y[-1], 2 * np.exp(0.1 * solution.t))
 
+    def test_model_solver_dae__with_external(self):
+        # Create model
+        model = pybamm.BaseModel()
+        domain = ["negative electrode", "separator", "positive electrode"]
+        var1 = pybamm.Variable("var1", domain=domain)
+        var2 = pybamm.Variable("var2", domain=domain)
+        model.rhs = {var1: -var2}
+        model.initial_conditions = {var1: 1}
+        model.external_variables = [var2]
+        model.variables = {"var1": var1, "var2": var2}
+        # No need to set parameters; can use base discretisation (no spatial
+        # operators)
+
+        # create discretisation
+        mesh = get_mesh_for_testing()
+        spatial_methods = {"macroscale": pybamm.FiniteVolume()}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+        disc.process_model(model)
+        # Solve
+        solver = pybamm.ScikitsDaeSolver(rtol=1e-8, atol=1e-8)
+        t_eval = np.linspace(0, 10, 100)
+        solution = solver.solve(model, t_eval, external_variables={"var2": 0.5})
+        np.testing.assert_allclose(solution.y[0], 1 - 0.5 * solution.t, rtol=1e-06)
+
     def test_solve_ode_model_with_dae_solver_casadi(self):
         model = pybamm.BaseModel()
         model.convert_to_format = "casadi"