From bbcfc22368ee9096ab04d660111b385588ae0567 Mon Sep 17 00:00:00 2001 From: Ferran Brosa Planella Date: Wed, 30 Dec 2020 18:06:54 +0100 Subject: [PATCH] #976 added error for extrapolation with CasADI --- pybamm/models/event.py | 1 + pybamm/parameters/parameter_values.py | 29 ++++++++++ pybamm/solvers/base_solver.py | 7 +++ pybamm/solvers/casadi_solver.py | 56 +++++++++++++++++++ tests/unit/test_solvers/test_casadi_solver.py | 22 +++++++- 5 files changed, 112 insertions(+), 3 deletions(-) diff --git a/pybamm/models/event.py b/pybamm/models/event.py index 578710e952..1fa86b8793 100644 --- a/pybamm/models/event.py +++ b/pybamm/models/event.py @@ -16,6 +16,7 @@ class EventType(Enum): TERMINATION = 0 DISCONTINUITY = 1 + INTERPOLANT_EXTRAPOLATION = 2 class Event: diff --git a/pybamm/parameters/parameter_values.py b/pybamm/parameters/parameter_values.py index 97eeb420c2..cac423ebf3 100644 --- a/pybamm/parameters/parameter_values.py +++ b/pybamm/parameters/parameter_values.py @@ -83,6 +83,7 @@ def __init__(self, values=None, chemistry=None): # Initialise empty _processed_symbols dict (for caching) self._processed_symbols = {} + self.parameter_events = [] def __getitem__(self, key): return self._dict_items[key] @@ -410,6 +411,17 @@ def process_model(self, unprocessed_model, inplace=True): event.name, self.process_symbol(event.expression), event.event_type ) ) + + for event in self.parameter_events: + pybamm.logger.debug( + "Processing parameters for event'{}''".format(event.name) + ) + new_events.append( + pybamm.Event( + event.name, self.process_symbol(event.expression), event.event_type + ) + ) + model.events = new_events # Set external variables @@ -545,6 +557,23 @@ def _process_symbol(self, symbol): # to create an Interpolant name, data = function_name function = pybamm.Interpolant(data, *new_children, name=name) + # Define event to catch extrapolation. In these events the sign is + # important: it should be positive inside of the range and negative + # outside of it + self.parameter_events.append( + pybamm.Event( + "Interpolant {} lower bound".format(name), + new_children[0] - min(function.x), + pybamm.EventType.INTERPOLANT_EXTRAPOLATION, + ) + ) + self.parameter_events.append( + pybamm.Event( + "Interpolant {} upper bound".format(name), + max(function.x) - new_children[0], + pybamm.EventType.INTERPOLANT_EXTRAPOLATION, + ) + ) elif isinstance(function_name, numbers.Number): # If the "function" is provided is actually a scalar, return a Scalar # object instead of throwing an error. diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py index 60533452f3..03c84d7a50 100644 --- a/pybamm/solvers/base_solver.py +++ b/pybamm/solvers/base_solver.py @@ -361,6 +361,12 @@ def report(string): if event.event_type == pybamm.EventType.TERMINATION ] + interpolant_extrapolation_events_eval = [ + process(event.expression, "event", use_jacobian=False)[1] + for event in model.events + if event.event_type == pybamm.EventType.INTERPOLANT_EXTRAPOLATION + ] + # discontinuity events are evaluated before the solver is called, so don't need # to process them discontinuity_events_eval = [ @@ -376,6 +382,7 @@ def report(string): model.jac_algebraic_eval = jac_algebraic model.terminate_events_eval = terminate_events_eval model.discontinuity_events_eval = discontinuity_events_eval + model.interpolant_extrapolation_events_eval = interpolant_extrapolation_events_eval # Calculate initial conditions model.y0 = init_eval(inputs) diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py index d05e6283f9..bb43dd9724 100644 --- a/pybamm/solvers/casadi_solver.py +++ b/pybamm/solvers/casadi_solver.py @@ -71,6 +71,7 @@ def __init__( root_tol=1e-6, max_step_decrease_count=5, dt_max=None, + extrap_tol=1e-3, extra_options_setup=None, extra_options_call=None, ): @@ -88,6 +89,7 @@ def __init__( self.extra_options_setup = extra_options_setup or {} self.extra_options_call = extra_options_call or {} + self.extrap_tol = extrap_tol self.name = "CasADi solver with '{}' mode".format(mode) @@ -143,6 +145,32 @@ def _integrate(self, model, t_eval, inputs=None): [event(t, y0, inputs) for event in model.terminate_events_eval] ) ) + + extrap_event = [ + event(t, y0, inputs) + for event in model.interpolant_extrapolation_events_eval + ] + + if extrap_event: + if (np.concatenate(extrap_event) < self.extrap_tol).any(): + extrap_event_names = [] + for event in model.events: + if ( + event.event_type + == pybamm.EventType.INTERPOLANT_EXTRAPOLATION + and ( + event.expression.evaluate(t, y0, inputs=inputs,) + < self.extrap_tol + ).any() + ): + extrap_event_names.append(event.name[12:]) + + raise pybamm.SolverError( + "CasADI solver failed because the following interpolation bounds were exceeded: {}".format( + extrap_event_names + ) + ) + pybamm.logger.info("Start solving {} with {}".format(model.name, self.name)) if self.mode == "safe without grid": @@ -217,6 +245,34 @@ def _integrate(self, model, t_eval, inputs=None): ] ) ) + + extrap_event = [ + event(t, current_step_sol.y[:, -1], inputs=inputs) + for event in model.interpolant_extrapolation_events_eval + ] + + if extrap_event: + if (np.concatenate(extrap_event) < self.extrap_tol).any(): + extrap_event_names = [] + for event in model.events: + if ( + event.event_type + == pybamm.EventType.INTERPOLANT_EXTRAPOLATION + and ( + event.expression.evaluate( + t, current_step_sol.y[:, -1], inputs=inputs + ) + < self.extrap_tol + ).any() + ): + extrap_event_names.append(event.name[12:]) + + raise pybamm.SolverError( + "CasADI solver failed because the following interpolation bounds were exceeded: {}".format( + extrap_event_names + ) + ) + # Exit loop if the sign of an event changes # Locate the event time using a root finding algorithm and # event state using interpolation. The solution is then truncated diff --git a/tests/unit/test_solvers/test_casadi_solver.py b/tests/unit/test_solvers/test_casadi_solver.py index 271d42be9d..2447b2a6be 100644 --- a/tests/unit/test_solvers/test_casadi_solver.py +++ b/tests/unit/test_solvers/test_casadi_solver.py @@ -426,9 +426,22 @@ def test_interpolant_extrapolate(self): model = pybamm.lithium_ion.DFN() param = pybamm.ParameterValues(chemistry=pybamm.parameter_sets.Chen2020) experiment = pybamm.Experiment( - ["Discharge at 1C until 2.5 V", "Rest for 2 hours",], period="5 seconds" + ["Charge at 1C until 4.6 V"], period="10 seconds" ) + param["Upper voltage cut-off [V]"] = 4.8 + + sim = pybamm.Simulation( + model, + parameter_values=param, + experiment=experiment, + solver=pybamm.CasadiSolver(mode="safe", dt_max=0.001, extra_options_setup={"max_num_steps": 500}), + ) + with self.assertRaisesRegex( + pybamm.SolverError, "interpolation bounds" + ): + sim.solve() + ci = param["Initial concentration in positive electrode [mol.m-3]"] param["Initial concentration in positive electrode [mol.m-3]"] = 0.8 * ci @@ -436,9 +449,12 @@ def test_interpolant_extrapolate(self): model, parameter_values=param, experiment=experiment, - solver=pybamm.CasadiSolver(mode="safe"), + solver=pybamm.CasadiSolver(mode="safe", dt_max=0.05), ) - sim.solve() + with self.assertRaisesRegex( + pybamm.SolverError, "interpolation bounds" + ): + sim.solve() class TestCasadiSolverSensitivity(unittest.TestCase):