Repository: pybamm-team/PyBaMM
Pull request #787 (2020-01-16 16:49:31)

Issue 775 casadi external
# Description

Fixes external variables with the casadi solver.
Will leave this issue open as there is more to do.

## Type of change

Please add a line in the relevant section of [CHANGELOG.md](https://github.com/pybamm-team/PyBaMM/blob/master/CHANGELOG.md) to document the change (include PR #) - note reverse order of PR #s. If necessary, also add to the list of breaking changes.

- [ ] New feature (non-breaking change which adds functionality)
- [ ] Optimization (back-end change that speeds up the code)
- [x] Bug fix (non-breaking change which fixes an issue)


# Key checklist:

- [x] No style issues: `$ flake8`
- [x] All tests pass: `$ python run-tests.py --unit`
- [x] The documentation builds: `$ cd docs` and then `$ make clean; make html`

You can run all three at once, using `$ python run-tests.py --quick`.

## Further checks:

- [x] Code is commented, particularly in hard-to-understand areas
- [x] Tests added that prove fix is effective or that feature works



Patch:
diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 7e93be8b5..cb465c100 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -49,7 +49,7 @@ class BaseSolver(object):
     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 @@ class BaseSolver(object):
             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 @@ class BaseSolver(object):
         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:
diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py
index 052c0d4ad..973d2b819 100644
--- a/pybamm/solvers/casadi_solver.py
+++ b/pybamm/solvers/casadi_solver.py
@@ -68,7 +68,7 @@ class CasadiSolver(pybamm.DaeSolver):
         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 @@ class CasadiSolver(pybamm.DaeSolver):
             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 @@ class CasadiSolver(pybamm.DaeSolver):
         """
         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 @@ class CasadiSolver(pybamm.DaeSolver):
                     # 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 @@ class CasadiSolver(pybamm.DaeSolver):
             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,14 +259,22 @@ class CasadiSolver(pybamm.DaeSolver):
         # set up and solve
         t = casadi.MX.sym("t")
         u = casadi.vertcat(*[x for x in inputs.values()])
-        y_diff = casadi.MX.sym("y_diff", rhs(0, y0, u).shape[0])
+        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:
-            problem.update({"ode": rhs(t, y_diff, u)})
+            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 = casadi.MX.sym("y_alg", algebraic(0, y0, u).shape[0])
-            y = casadi.vertcat(y_diff, y_alg)
-            problem.update({"z": y_alg, "ode": rhs(t, y, u), "alg": algebraic(t, y, u)})
+            y_alg = casadi.MX.sym("y_alg", algebraic(0, y0_w_ext, u).shape[0])
+            y_casadi_w_ext = casadi.vertcat(y_diff, y_alg, y_ext[len(y0) :])
+            problem.update(
+                {
+                    "z": y_alg,
+                    "ode": rhs(t, y_casadi_w_ext, u),
+                    "alg": algebraic(t, y_casadi_w_ext, u),
+                }
+            )
         integrator = casadi.integrator("F", self.method, problem, options)
         try:
             # Try solving
diff --git a/pybamm/solvers/ode_solver.py b/pybamm/solvers/ode_solver.py
index 135d93dfb..ffd737a2c 100644
--- a/pybamm/solvers/ode_solver.py
+++ b/pybamm/solvers/ode_solver.py
@@ -188,11 +188,8 @@ class OdeSolver(pybamm.BaseSolver):
         inputs = inputs or {}
         u_casadi = {name: casadi.MX.sym(name) for name in inputs.keys()}
 
-        if self.y_pad is not None:
-            y_ext = casadi.MX.sym("y_ext", len(self.y_pad))
-            y_casadi_w_ext = casadi.vertcat(y_casadi, y_ext)
-        else:
-            y_casadi_w_ext = y_casadi
+        y_ext = casadi.MX.sym("y_ext", len(self.y_pad))
+        y_casadi_w_ext = casadi.vertcat(y_casadi, y_ext)
 
         pybamm.logger.info("Converting RHS to CasADi")
         concatenated_rhs = model.concatenated_rhs.to_casadi(


Test patch:
diff --git a/tests/unit/test_solvers/test_casadi_solver.py b/tests/unit/test_solvers/test_casadi_solver.py
index 01b5b2f9e..cda4d6140 100644
--- a/tests/unit/test_solvers/test_casadi_solver.py
+++ b/tests/unit/test_solvers/test_casadi_solver.py
@@ -110,7 +110,7 @@ class TestCasadiSolver(unittest.TestCase):
         disc = pybamm.Discretisation(mesh, spatial_methods)
         disc.process_model(model)
         # Solve
-        solver = pybamm.CasadiSolver(rtol=1e-8, atol=1e-8, method="idas")
+        solver = pybamm.CasadiSolver(mode="fast", rtol=1e-8, atol=1e-8, method="idas")
         t_eval = np.linspace(0, 1, 100)
         solution = solver.solve(model, t_eval)
         np.testing.assert_array_equal(solution.t, t_eval)
@@ -214,6 +214,30 @@ class TestCasadiSolver(unittest.TestCase):
         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")
diff --git a/tests/unit/test_solvers/test_ode_solver.py b/tests/unit/test_solvers/test_ode_solver.py
index 5469e4f72..c079681c2 100644
--- a/tests/unit/test_solvers/test_ode_solver.py
+++ b/tests/unit/test_solvers/test_ode_solver.py
@@ -14,9 +14,13 @@ class TestOdeSolver(unittest.TestCase):
     def test_wrong_solver(self):
         # Create model
         model = pybamm.BaseModel()
-        var = pybamm.Variable("var")
-        model.rhs = {var: var}
-        model.algebraic = {var: var - 1}
+        var1 = pybamm.Variable("var1")
+        var2 = pybamm.Variable("var2")
+        model.rhs = {var1: var1}
+        model.algebraic = {var2: var2 - 1}
+        model.initial_conditions = {var1: 0, var2: 1}
+        disc = pybamm.Discretisation()
+        disc.process_model(model)
 
         # test errors
         solver = pybamm.OdeSolver()
diff --git a/tests/unit/test_solvers/test_scikits_solvers.py b/tests/unit/test_solvers/test_scikits_solvers.py
index c2259b6a7..3d512cfdf 100644
--- a/tests/unit/test_solvers/test_scikits_solvers.py
+++ b/tests/unit/test_solvers/test_scikits_solvers.py
@@ -745,6 +745,30 @@ class TestScikitsSolvers(unittest.TestCase):
             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"
diff --git a/tests/unit/test_solvers/test_scipy_solver.py b/tests/unit/test_solvers/test_scipy_solver.py
index ff10eadc8..2dde32f93 100644
--- a/tests/unit/test_solvers/test_scipy_solver.py
+++ b/tests/unit/test_solvers/test_scipy_solver.py
@@ -297,6 +297,30 @@ class TestScipySolver(unittest.TestCase):
         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))
 
+    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.ScipySolver(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_event_with_casadi(self):
         # Create model
         model = pybamm.BaseModel()