From e71a4eef9adf77ae93cb8ca483c99d26cb6d5044 Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Tue, 16 Mar 2021 14:34:46 -0400 Subject: [PATCH] #1082 basic 'fast with events' algorithm working but needs fine tuning --- examples/scripts/DFN.py | 22 +- .../scripts/experimental_protocols/cccv.py | 30 +- .../graphite_ocp_PeymanMPM.py | 11 +- .../NMC_UMBL_Mohtat2020/NMC_ocp_PeymanMPM.py | 9 +- .../full_battery_models/base_battery_model.py | 14 +- .../submodels/interface/base_interface.py | 5 + pybamm/parameters/lithium_ion_parameters.py | 14 +- pybamm/simulation.py | 15 +- pybamm/solvers/base_solver.py | 14 +- pybamm/solvers/casadi_solver.py | 328 ++++++++++-------- 10 files changed, 264 insertions(+), 198 deletions(-) diff --git a/examples/scripts/DFN.py b/examples/scripts/DFN.py index 33fffacef5..bc0b26fd11 100644 --- a/examples/scripts/DFN.py +++ b/examples/scripts/DFN.py @@ -8,12 +8,14 @@ pybamm.set_logging_level("INFO") # load model -model = pybamm.lithium_ion.DFN() +model = pybamm.lithium_ion.DFN() # {"operating mode": "power"}) # create geometry geometry = model.default_geometry # load parameter values and process model and geometry param = model.default_parameter_values +# param.update({"Power function [W]": 3.5}, check_already_exists=False) +param["Current function [A]"] /= 10 param.process_geometry(geometry) param.process_model(model) @@ -27,21 +29,23 @@ disc.process_model(model) # solve model -t_eval = np.linspace(0, 5000, 100) -solver = pybamm.CasadiSolver(mode="fast", atol=1e-6, rtol=1e-3) -solution = solver.solve(model, t_eval) +t_eval = np.linspace(0, 5000 * 10, 100) +solver = pybamm.CasadiSolver(mode="safe", atol=1e-6, rtol=1e-6) +solution1 = solver.solve(model, t_eval) +solver = pybamm.CasadiSolver(mode="fast with events", atol=1e-6, rtol=1e-6) +solution2 = solver.solve(model, t_eval) # plot plot = pybamm.QuickPlot( - solution, + [solution1, solution2], [ - "Negative particle concentration [mol.m-3]", + # "Negative particle concentration [mol.m-3]", "Electrolyte concentration [mol.m-3]", - "Positive particle concentration [mol.m-3]", + # "Positive particle concentration [mol.m-3]", "Current [A]", "Negative electrode potential [V]", - "Electrolyte potential [V]", - "Positive electrode potential [V]", + # "Electrolyte potential [V]", + # "Positive electrode potential [V]", "Terminal voltage [V]", ], time_unit="seconds", diff --git a/examples/scripts/experimental_protocols/cccv.py b/examples/scripts/experimental_protocols/cccv.py index 7d824c6c94..8a75b77ebd 100644 --- a/examples/scripts/experimental_protocols/cccv.py +++ b/examples/scripts/experimental_protocols/cccv.py @@ -4,7 +4,7 @@ import pybamm import matplotlib.pyplot as plt -pybamm.set_logging_level("VERBOSE") +pybamm.set_logging_level("INFO") experiment = pybamm.Experiment( [ ( @@ -19,24 +19,24 @@ ) model = pybamm.lithium_ion.SPM() sim = pybamm.Simulation( - model, experiment=experiment, solver=pybamm.CasadiSolver(dt_max=1e4) + model, experiment=experiment, solver=pybamm.CasadiSolver("fast with events") ) sim.solve() # Plot voltages from the discharge segments only -fig, ax = plt.subplots() -for i in range(3): - # Extract sub solutions - sol = sim.solution.cycles[i] - # Extract variables - t = sol["Time [h]"].entries - V = sol["Terminal voltage [V]"].entries - # Plot - ax.plot(t - t[0], V, label="Discharge {}".format(i + 1)) - ax.set_xlabel("Time [h]") - ax.set_ylabel("Voltage [V]") - ax.set_xlim([0, 10]) -ax.legend(loc="lower left") +# fig, ax = plt.subplots() +# for i in range(3): +# # Extract sub solutions +# sol = sim.solution.cycles[i] +# # Extract variables +# t = sol["Time [h]"].entries +# V = sol["Terminal voltage [V]"].entries +# # Plot +# ax.plot(t - t[0], V, label="Discharge {}".format(i + 1)) +# ax.set_xlabel("Time [h]") +# ax.set_ylabel("Voltage [V]") +# ax.set_xlim([0, 10]) +# ax.legend(loc="lower left") # Save time, voltage, current, discharge capacity, temperature, and electrolyte # concentration to csv and matlab formats diff --git a/pybamm/input/parameters/lithium-ion/negative_electrodes/graphite_UMBL_Mohtat2020/graphite_ocp_PeymanMPM.py b/pybamm/input/parameters/lithium-ion/negative_electrodes/graphite_UMBL_Mohtat2020/graphite_ocp_PeymanMPM.py index 02a5bba688..4bfe02dd87 100644 --- a/pybamm/input/parameters/lithium-ion/negative_electrodes/graphite_UMBL_Mohtat2020/graphite_ocp_PeymanMPM.py +++ b/pybamm/input/parameters/lithium-ion/negative_electrodes/graphite_UMBL_Mohtat2020/graphite_ocp_PeymanMPM.py @@ -25,10 +25,7 @@ def graphite_ocp_PeymanMPM(sto): return u_eq -# if __name__ == "__main__": # pragma: no cover -# import matplotlib.pyplot as plt -# import numpy as np - -# x = np.linspace(0, 1) -# plt.plot(x, graphite_ocp_PeymanMPM(x)) -# plt.show() +if __name__ == "__main__": # pragma: no cover + x = pybamm.linspace(1e-10, 1 - 1e-10, 1000) + # pybamm.plot(x, graphite_ocp_PeymanMPM(x)) + pybamm.plot(x, -1e-8 * pybamm.log(x / (1 - x))) diff --git a/pybamm/input/parameters/lithium-ion/positive_electrodes/NMC_UMBL_Mohtat2020/NMC_ocp_PeymanMPM.py b/pybamm/input/parameters/lithium-ion/positive_electrodes/NMC_UMBL_Mohtat2020/NMC_ocp_PeymanMPM.py index 8cde5060df..b2eb8f54fd 100644 --- a/pybamm/input/parameters/lithium-ion/positive_electrodes/NMC_UMBL_Mohtat2020/NMC_ocp_PeymanMPM.py +++ b/pybamm/input/parameters/lithium-ion/positive_electrodes/NMC_UMBL_Mohtat2020/NMC_ocp_PeymanMPM.py @@ -30,9 +30,6 @@ def NMC_ocp_PeymanMPM(sto): return u_eq -# if __name__ == "__main__": # pragma: no cover -# import matplotlib.pyplot as plt - -# x = np.linspace(0, 1) -# plt.plot(x, NMC_ocp_PeymanMPM(x)) -# plt.show() +# if __name__ == "__main__": # pragma: no cover +# x = pybamm.linspace(0, 1) +# pybamm.plot(x, NMC_ocp_PeymanMPM(x)) diff --git a/pybamm/models/full_battery_models/base_battery_model.py b/pybamm/models/full_battery_models/base_battery_model.py index 83942c997d..bb206670c4 100644 --- a/pybamm/models/full_battery_models/base_battery_model.py +++ b/pybamm/models/full_battery_models/base_battery_model.py @@ -899,13 +899,13 @@ def set_voltage_variables(self): ) # Cut-off voltage - # self.events.append( - # pybamm.Event( - # "Minimum voltage", - # V - self.param.voltage_low_cut, - # pybamm.EventType.TERMINATION, - # ) - # ) + self.events.append( + pybamm.Event( + "Minimum voltage", + V - self.param.voltage_low_cut, + pybamm.EventType.TERMINATION, + ) + ) self.events.append( pybamm.Event( "Maximum voltage", diff --git a/pybamm/models/submodels/interface/base_interface.py b/pybamm/models/submodels/interface/base_interface.py index e80ce1e1ed..fd9cbef588 100644 --- a/pybamm/models/submodels/interface/base_interface.py +++ b/pybamm/models/submodels/interface/base_interface.py @@ -72,6 +72,11 @@ def _get_exchange_current_density(self, variables): c_s_surf = c_s_surf.orphans[0] c_e = c_e.orphans[0] T = T.orphans[0] + + tol = 1e-8 + c_e = pybamm.maximum(tol, c_e) + c_s_surf = pybamm.maximum(tol, pybamm.minimum(c_s_surf, 1 - tol)) + if self.domain == "Negative": j0 = self.param.j0_n(c_e, c_s_surf, T) / self.param.C_r_n elif self.domain == "Positive": diff --git a/pybamm/parameters/lithium_ion_parameters.py b/pybamm/parameters/lithium_ion_parameters.py index 5400c37c14..0720d14812 100644 --- a/pybamm/parameters/lithium_ion_parameters.py +++ b/pybamm/parameters/lithium_ion_parameters.py @@ -367,12 +367,20 @@ def U_n_dimensional(self, sto, T): """Dimensional open-circuit potential in the negative electrode [V]""" inputs = {"Negative particle stoichiometry": sto} u_ref = pybamm.FunctionParameter("Negative electrode OCP [V]", inputs) + # add a term to ensure that the OCP goes to infinity at 0 and -infinity at 1 + # this will not affect the OCP for most values of sto + # see #1435 + u_ref -= 1e-6 * pybamm.log(sto / (1 - sto)) return u_ref + (T - self.T_ref) * self.dUdT_n_dimensional(sto) def U_p_dimensional(self, sto, T): """Dimensional open-circuit potential in the positive electrode [V]""" inputs = {"Positive particle stoichiometry": sto} u_ref = pybamm.FunctionParameter("Positive electrode OCP [V]", inputs) + # add a term to ensure that the OCP goes to infinity at 0 and -infinity at 1 + # this will not affect the OCP for most values of sto + # see #1435 + u_ref -= 1e-6 * pybamm.log(sto / (1 - sto)) return u_ref + (T - self.T_ref) * self.dUdT_p_dimensional(sto) def dUdT_n_dimensional(self, sto): @@ -880,13 +888,11 @@ def R_p(self, x): return self.R_p_dimensional(x_dim) / self.R_p_typ def c_n_init(self, x): - """Dimensionless initial concentration as a function of dimensionless position x - """ + """Dimensionless initial concentration as a function of dimensionless position x""" return self.c_n_init_dimensional(x) / self.c_n_max def c_p_init(self, x): - """Dimensionless initial concentration as a function of dimensionless position x - """ + """Dimensionless initial concentration as a function of dimensionless position x""" return self.c_p_init_dimensional(x) / self.c_p_max def rho(self, T): diff --git a/pybamm/simulation.py b/pybamm/simulation.py index e5c4f0e3a3..38fca7957b 100644 --- a/pybamm/simulation.py +++ b/pybamm/simulation.py @@ -172,12 +172,14 @@ def set_up_experiment(self, model, experiment): for op, events in zip(experiment.operating_conditions, experiment.events): if op[1] in ["A", "C"]: # Update inputs for constant current + capacity = self._parameter_values["Nominal cell capacity [A.h]"] if op[1] == "A": I = op[0] + Crate = I / capacity else: # Scale C-rate with capacity to obtain current - capacity = self._parameter_values["Nominal cell capacity [A.h]"] - I = op[0] * capacity + Crate = op[0] + I = Crate * capacity operating_inputs = { "Current switch": 1, "Voltage switch": 0, @@ -238,8 +240,13 @@ def set_up_experiment(self, model, experiment): # Add time to the experiment times dt = op[2] if dt is None: - # max simulation time: 1 week - dt = 7 * 24 * 3600 + if op[1] in ["A", "C"]: + # Current control: max simulation time: 3 * max simulation time + # based on C-rate + dt = 3 / abs(Crate) + else: + # max simulation time: 1 week + dt = 7 * 24 * 3600 self._experiment_times.append(dt) # Set up model for experiment diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py index d1d31ee4df..12bc66ddf4 100644 --- a/pybamm/solvers/base_solver.py +++ b/pybamm/solvers/base_solver.py @@ -378,15 +378,17 @@ def report(string): # events' mode of the casadi solver # see #1082 k = 20 - if ( - event.name in ["Minimum voltage", "Maximum voltage"] - or "[experiment]" in event.name - ): + if "voltage" in event.name.lower(): init_sign = float( np.sign(event_eval(0, model.y0, inputs_stacked)) ) - event_sigmoid = pybamm.sigmoid( - 0, init_sign * event.expression, k + # We create a sigmoid for each event which will multiply the + # rhs. Doing * 2 - 1 ensures that when the event is crossed, + # the sigmoid is zero. Hence the rhs is zero and the solution + # stays constant for the rest of the simulation period + # We can then cut off the part after the event was crossed + event_sigmoid = ( + pybamm.sigmoid(0, init_sign * event.expression, k) * 2 - 1 ) event_casadi = process( event_sigmoid, "event", use_jacobian=False diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py index 7d344be997..312a21686c 100644 --- a/pybamm/solvers/casadi_solver.py +++ b/pybamm/solvers/casadi_solver.py @@ -21,6 +21,8 @@ class CasadiSolver(pybamm.BaseSolver): - "fast": perform direct integration, without accounting for events. \ Recommended when simulating a drive cycle or other simulation where \ no events should be triggered. + - "fast with events": perform direct integration of the whole timespan, \ + then go back and check where events were crossed. Experimental only. - "safe": perform step-and-check integration in global steps of size \ dt_max, checking whether events have been triggered. Recommended for \ simulations of a full charge or discharge. @@ -125,6 +127,23 @@ def _integrate(self, model, t_eval, inputs_dict=None): # convert inputs to casadi format inputs = casadi.vertcat(*[x for x in inputs_dict.values()]) + # Calculate initial event signs needed for some of the modes + if ( + has_symbolic_inputs is False + and self.mode != "fast" + and model.terminate_events_eval + ): + init_event_signs = np.sign( + np.concatenate( + [ + event(t_eval[0], model.y0, inputs) + for event in model.terminate_events_eval + ] + ) + ) + else: + init_event_signs = np.sign([]) + if has_symbolic_inputs: # Create integrator without grid to avoid having to create several times self.create_integrator(model, inputs) @@ -133,38 +152,31 @@ def _integrate(self, model, t_eval, inputs_dict=None): ) solution.termination = "final time" return solution - elif self.mode == "fast" or not model.events: + elif self.mode in ["fast", "fast with events"] or not model.events: if not model.events: pybamm.logger.info("No events found, running fast mode") + if self.mode == "fast with events": + # Create the integrator with an event switch that will set the rhs to + # zero when voltage limits are crossed + use_event_switch = True + else: + use_event_switch = False # Create an integrator with the grid (we just need to do this once) - self.create_integrator(model, inputs, t_eval) - solution = self._run_integrator( - model, model.y0, inputs_dict, inputs, t_eval + self.create_integrator( + model, inputs, t_eval, use_event_switch=use_event_switch ) - solution.termination = "final time" - return solution - elif self.mode == "fast with events": - # This is not working but keeping in case we can find a way to make it - # work - self.create_integrator(model, inputs, t_eval, use_event_switch=True) solution = self._run_integrator( model, model.y0, inputs_dict, inputs, t_eval ) - solution.termination = "final time" + # Check if the sign of an event changes, if so find an accurate + # termination point and exit + solution = self._solve_for_event(solution, init_event_signs) return solution elif self.mode in ["safe", "safe without grid"]: y0 = model.y0 # Step-and-check t = t_eval[0] t_f = t_eval[-1] - if model.terminate_events_eval: - init_event_signs = np.sign( - np.concatenate( - [event(t, y0, inputs) for event in model.terminate_events_eval] - ) - ) - else: - init_event_signs = np.sign([]) pybamm.logger.debug( "Start solving {} with {}".format(model.name, self.name) @@ -233,135 +245,171 @@ def _integrate(self, model, t_eval, inputs_dict=None): t * model.timescale_eval, dt_max * model.timescale_eval ) ) - # Check most recent y to see if any events have been crossed - if model.terminate_events_eval: - new_event_signs = np.sign( - np.concatenate( - [ - event(t, current_step_sol.all_ys[-1][:, -1], inputs) - for event in model.terminate_events_eval - ] - ) - ) + # Check if the sign of an event changes, if so find an accurate + # termination point and exit + current_step_sol = self._solve_for_event( + current_step_sol, init_event_signs + ) + # assign temporary solve time + current_step_sol.solve_time = np.nan + # append solution from the current step to solution + solution = solution + current_step_sol + if current_step_sol.termination == "event": + break else: - new_event_signs = np.sign([]) + # update time + t = t_window[-1] + # update y0 + y0 = solution.all_ys[-1][:, -1] + return solution - if model.interpolant_extrapolation_events_eval: - extrap_event = [ - event(t, current_step_sol.y[:, -1], inputs=inputs) - for event in model.interpolant_extrapolation_events_eval - ] + def _solve_for_event(self, coarse_solution, init_event_signs): + """ + Check if the sign of an event changes, if so find an accurate + termination point and exit - 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].full(), - 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: {}. You may need " - "to provide additional interpolation points outside " - "these bounds.".format(extrap_event_names) - ) + Locate the event time using a root finding algorithm and + event state using interpolation. The solution is then truncated + so that only the times up to the event are returned + """ + model = coarse_solution.all_models[-1] + inputs_dict = coarse_solution.all_inputs[-1] + inputs = casadi.vertcat(*[x for x in inputs_dict.values()]) - # 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 - # so that only the times up to the event are returned - if (new_event_signs != init_event_signs).any(): - # get the index of the events that have been crossed - event_ind = np.where(new_event_signs != init_event_signs)[0] - active_events = [model.terminate_events_eval[i] for i in event_ind] - - # solve again with a more dense t_window - if len(t_window) < 10: - t_window_dense = np.linspace(t_window[0], t_window[-1], 10) - if self.mode == "safe": - self.create_integrator(model, inputs, t_window_dense) - current_step_sol = self._run_integrator( - model, y0, inputs_dict, inputs, t_window_dense - ) - integration_time = current_step_sol.integration_time + # Check for interpolant extrapolations + if model.interpolant_extrapolation_events_eval: + extrap_event = [ + event(t, coarse_solution.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, + coarse_solution.y[:, -1].full(), + inputs=inputs, + ) + < self.extrap_tol + ).any() + ): + extrap_event_names.append(event.name[12:]) - # create interpolant to evaluate y in the current integration - # window - y_sol = interp1d( - current_step_sol.t, current_step_sol.y, kind="cubic" + raise pybamm.SolverError( + "CasADI solver failed because the following " + "interpolation bounds were exceeded: {}. You may need " + "to provide additional interpolation points outside " + "these bounds.".format(extrap_event_names) ) - # loop over events to compute the time at which they were triggered - t_events = [None] * len(active_events) - for i, event in enumerate(active_events): - - def event_fun(t): - return event(t, y_sol(t), inputs) - - if np.isnan(event_fun(current_step_sol.t[-1])[0]): - # bracketed search fails if f(a) or f(b) is NaN, so we - # need to find the times for which we can evaluate the event - times = [ - t - for t in current_step_sol.t - if event_fun(t)[0] == event_fun(t)[0] - ] - else: - times = current_step_sol.t - # skip if sign hasn't changed - if np.sign(event_fun(times[0])) != np.sign( - event_fun(times[-1]) - ): - t_events[i] = brentq( - lambda t: event_fun(t), times[0], times[-1] - ) - else: - t_events[i] = np.nan - - # t_event is the earliest event triggered - t_event = np.nanmin(t_events) - y_event = y_sol(t_event) + def find_t_event(sol, interpolant_kind): - # call interpolant at times in t_eval and create solution - t_window = np.concatenate( - ([t], t_eval[(t_eval > t) & (t_eval < t_event)]) - ) - y_sol_casadi = casadi.DM(y_sol(t_window)) - current_step_sol = pybamm.Solution( - t_window, y_sol_casadi, model, inputs_dict + # Check most recent y to see if any events have been crossed + if model.terminate_events_eval: + y_last = sol.all_ys[-1][:, -1] + new_event_signs = np.sign( + np.concatenate( + [ + event(sol.t[-1], y_last, inputs) + for event in model.terminate_events_eval + ] ) - current_step_sol.integration_time = integration_time - # assign temporary solve time - current_step_sol.solve_time = np.nan + ) + else: + new_event_signs = np.sign([]) + + # Return None if no events have been triggered + if (new_event_signs == init_event_signs).all(): + return None, None + + # get the index of the events that have been crossed + event_ind = np.where(new_event_signs != init_event_signs)[0] + active_events = [model.terminate_events_eval[i] for i in event_ind] + + # create interpolant to evaluate y in the current integration + # window + y_sol = interp1d(sol.t, sol.y, kind=interpolant_kind) + + # loop over events to compute the time at which they were triggered + t_events = [None] * len(active_events) + for i, event in enumerate(active_events): + init_event_sign = init_event_signs[event_ind[i]][0] + + def event_fun(t): + # We take away 1e-5 to deal with the case where the event sits + # exactly on zero, as can happen when the event switch is used + # (fast with events mode) + return init_event_sign * event(t, y_sol(t), inputs) - 1e-5 + + if np.isnan(event_fun(sol.t[-1])[0]): + # bracketed search fails if f(a) or f(b) is NaN, so we + # need to find the times for which we can evaluate the event + times = [t for t in sol.t if event_fun(t)[0] == event_fun(t)[0]] + else: + times = sol.t + # skip if sign hasn't changed + if np.sign(event_fun(times[0])) != np.sign(event_fun(times[-1])): + t_events[i] = brentq(lambda t: event_fun(t), times[0], times[-1]) + else: + t_events[i] = np.nan - # append solution from the current step to solution - solution = solution + current_step_sol - solution.termination = "event" - solution.t_event = np.array([t_event]) - solution.y_event = y_event[:, np.newaxis] + # t_event is the earliest event triggered + t_event = np.nanmin(t_events) + y_event = y_sol(t_event) - break - else: - # assign temporary solve time - current_step_sol.solve_time = np.nan - # append solution from the current step to solution - solution = solution + current_step_sol - # update time - t = t_window[-1] - # update y0 - y0 = solution.all_ys[-1][:, -1] - return solution + return t_event, y_event + + # Find the interval in which the event was triggered + t_event_coarse, _ = find_t_event(coarse_solution, "linear") + + # Return the existing solution if no events have been triggered + if t_event_coarse is None: + # Flag "final time" for termination + coarse_solution.termination = "final time" + return coarse_solution + + # If events have been triggered, we solve for a dense window in the interval + # where the event was triggered, then find the precise location of the event + event_idx = np.where(coarse_solution.t > t_event_coarse)[0][0] - 1 + t_window_event = coarse_solution.t[event_idx : event_idx + 2] + + # Solve again with a more dense t_window, starting from the start of the + # window where the event was triggered + t_window_event_dense = np.linspace(t_window_event[-2], t_window_event[-1], 100) + if self.mode in ["safe", "fast with events"]: + self.create_integrator(model, inputs, t_window_event_dense) + + y0 = coarse_solution.y[:, event_idx] + dense_step_sol = self._run_integrator( + model, y0, inputs_dict, inputs, t_window_event_dense + ) + + # Find the exact time at which the event was triggered + t_event, y_event = find_t_event(dense_step_sol, "cubic") + + # Return solution truncated at the first coarse event time + # Also assign t_event + t_sol = coarse_solution.t[: event_idx + 1] + y_sol = coarse_solution.y[:, : event_idx + 1] + solution = pybamm.Solution( + t_sol, + y_sol, + model, + inputs_dict, + np.array([t_event]), + y_event[:, np.newaxis], + "event", + ) + solution.integration_time = coarse_solution.integration_time + + # Flag "True" for termination + return solution def create_integrator(self, model, inputs, t_eval=None, use_event_switch=False): """ @@ -432,9 +480,10 @@ def create_integrator(self, model, inputs, t_eval=None, use_event_switch=False): p_with_tlims = p # define the event switch as the point when an event is crossed + # we don't do this for ODE models # see #1082 event_switch = 1 - if use_event_switch is True: + if use_event_switch is True and not algebraic(0, y0, p).is_empty(): for event in model.casadi_terminate_events: event_switch *= event(t_scaled, y_full, p) @@ -452,8 +501,7 @@ def create_integrator(self, model, inputs, t_eval=None, use_event_switch=False): problem.update( { "z": y_alg, - "alg": algebraic(t_scaled, y_full, p) * event_switch - + y_alg * (1 - event_switch), + "alg": algebraic(t_scaled, y_full, p), } ) integrator = casadi.integrator("F", method, problem, options)