Skip to content

Commit

Permalink
#1082 reformat casadi solver with explicit 'safe without grid' mode
Browse files Browse the repository at this point in the history
  • Loading branch information
@valentinsulzer
valentinsulzer committed Jul 20, 2020
1 parent bb098dd commit d24d4aa
Show file tree
Hide file tree
Showing 2 changed files with 61 additions and 40 deletions.
92 changes: 53 additions & 39 deletions pybamm/solvers/casadi_solver.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,9 @@ class CasadiSolver(pybamm.BaseSolver):
- "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.
- "safe without grid": perform step-and-check integration step-by-step. \
Takes more steps than "safe" mode, but doesn't require creating the grid \
each time, so may be faster. Experimental only.
rtol : float, optional
The relative tolerance for the solver (default is 1e-6).
atol : float, optional
Expand Down Expand Up @@ -71,12 +74,13 @@ def __init__(
extra_options_call=None,
):
super().__init__("problem dependent", rtol, atol, root_method, root_tol)
if mode in ["safe", "fast"]:
if mode in ["safe", "fast", "safe without grid"]:
self.mode = mode
else:
raise ValueError(
"invalid mode '{}'. Must be 'safe', for solving with events, "
"or 'fast', for solving quickly without events".format(mode)
"'fast', for solving quickly without events, or 'safe without grid' "
"(experimental)".format(mode)
)
self.max_step_decrease_count = max_step_decrease_count
self.dt_max = dt_max
Expand Down Expand Up @@ -113,22 +117,20 @@ def _integrate(self, model, t_eval, inputs=None):
inputs = casadi.vertcat(*[x for x in inputs.values()])

if has_symbolic_inputs:
# Create integrax`tor without grid to avoid having to create several times
self.get_integrator(model, inputs)
# Create integrator without grid to avoid having to create several times
self.create_integrator(model, inputs)
solution = self._run_integrator(model, model.y0, inputs, t_eval)
solution.termination = "final time"
return solution
elif self.mode == "fast" or not model.events:
if not model.events:
pybamm.logger.info("No events found, running fast mode")
# Create an integrator with the grid (we just need to do this once)
self.get_integrator(model, inputs, t_eval)
self.create_integrator(model, inputs, t_eval)
solution = self._run_integrator(model, model.y0, inputs, t_eval)
solution.termination = "final time"
return solution
elif self.mode == "safe":
# Create integrator without grid to avoid having to create several times
self.get_integrator(model, inputs)
elif self.mode in ["safe", "safe without grid"]:
y0 = model.y0
if isinstance(y0, casadi.DM):
y0 = y0.full().flatten()
Expand All @@ -142,9 +144,16 @@ def _integrate(self, model, t_eval, inputs=None):
)
pybamm.logger.info("Start solving {} with {}".format(model.name, self.name))

# Initialize solution
solution = pybamm.Solution(np.array([t]), y0[:, np.newaxis])
solution.solve_time = 0
if self.mode == "safe without grid":
# in "safe without grid" mode,
# create integrator once, without grid,
# to avoid having to create several times
self.create_integrator(model, inputs)
# Initialize solution
solution = pybamm.Solution(np.array([t]), y0[:, np.newaxis])
solution.solve_time = 0
else:
solution = None

# Try to integrate in global steps of size dt_max. Note: dt_max must
# be at least as big as the the biggest step in t_eval (multiplied
Expand Down Expand Up @@ -173,7 +182,9 @@ def _integrate(self, model, t_eval, inputs=None):
if len(t_window) == 1:
t_window = np.array([t, t + dt])

# integrator = self.get_integrator(model, t_window, inputs)
if self.mode == "safe":
# update integrator with the grid
self.create_integrator(model, inputs, t_window)
# Try to solve with the current global step, if it fails then
# halve the step size and try again.
try:
Expand Down Expand Up @@ -256,13 +267,17 @@ def event_fun(t):
if len(t_window) == 1:
t_window = np.array([t, t_event])

# integrator = self.get_integrator(model, t_window, inputs)
if self.mode == "safe":
self.create_integrator(model, inputs, t_window)
current_step_sol = self._run_integrator(model, y0, inputs, t_window)

# assign temporary solve time
current_step_sol.solve_time = np.nan
# append solution from the current step to solution
solution.append(current_step_sol)
if solution is None:
solution = current_step_sol
else:
# append solution from the current step to solution
solution.append(current_step_sol)
solution.termination = "event"
solution.t_event = t_event
solution.y_event = y_event
Expand All @@ -271,15 +286,18 @@ def event_fun(t):
else:
# assign temporary solve time
current_step_sol.solve_time = np.nan
# append solution from the current step to solution
solution.append(current_step_sol)
if solution is None:
solution = current_step_sol
else:
# append solution from the current step to solution
solution.append(current_step_sol)
# update time
t = t_window[-1]
# update y0
y0 = solution.y[:, -1]
return solution

def get_integrator(self, model, inputs, t_eval=None):
def create_integrator(self, model, inputs, t_eval=None):
"""
Method to create a casadi integrator object.
If t_eval is provided, the integrator uses t_eval to make the grid.
Expand All @@ -292,11 +310,10 @@ def get_integrator(self, model, inputs, t_eval=None):
# If we're not using the grid, we don't need to change the integrator
if use_grid is False:
return self.integrators[model]
# Otherwise, create new integrator with an updated (scaled) grid
# Otherwise, create new integrator with an updated grid
else:
method, problem, options = self.integrator_specs[model]
t_eval_scaled = (t_eval - t_eval[0]) / (t_eval[-1] - t_eval[0])
options["grid"] = t_eval_scaled
options["grid"] = t_eval
integrator = casadi.integrator("F", method, problem, options)
self.integrators[model] = (integrator, use_grid)
return integrator
Expand Down Expand Up @@ -326,17 +343,20 @@ def get_integrator(self, model, inputs, t_eval=None):
p = casadi.MX.sym("p", inputs.shape[0])
y_diff = casadi.MX.sym("y_diff", rhs(0, y0, p).shape[0])

# rescale time
t_min = casadi.MX.sym("t_min")
t_max = casadi.MX.sym("t_max")
t_scaled = t_min + (t_max - t_min) * t
# add time limits as inputs
p_with_tlims = casadi.vertcat(p, t_min, t_max)

# save (scaled) grid
if use_grid is True:
t_eval_scaled = (t_eval - t_eval[0]) / (t_eval[-1] - t_eval[0])
options.update({"grid": t_eval_scaled, "output_t0": True})
if use_grid is False:
# rescale time
t_min = casadi.MX.sym("t_min")
t_max = casadi.MX.sym("t_max")
t_scaled = t_min + (t_max - t_min) * t
# add time limits as inputs
p_with_tlims = casadi.vertcat(p, t_min, t_max)
else:
options.update({"grid": t_eval, "output_t0": True})
# Set dummy parameters for consistency with rescaled time
t_max = 1
t_min = 0
t_scaled = t
p_with_tlims = p

problem = {"t": t, "x": y_diff, "p": p_with_tlims}
if algebraic(0, y0, p).is_empty():
Expand Down Expand Up @@ -369,15 +389,9 @@ def _run_integrator(self, model, y0, inputs, t_eval):
try:
# Try solving
if use_grid is True:
t_min = t_eval[0]
t_max = t_eval[-1]
inputs_with_tlims = casadi.vertcat(inputs, t_min, t_max)
# Call the integrator once, with the grid
sol = integrator(
x0=y0_diff,
z0=y0_alg,
p=inputs_with_tlims,
**self.extra_options_call
x0=y0_diff, z0=y0_alg, p=inputs, **self.extra_options_call
)
y_sol = np.concatenate([sol["xf"].full(), sol["zf"].full()])
return pybamm.Solution(t_eval, y_sol)
Expand Down
9 changes: 8 additions & 1 deletion tests/unit/test_solvers/test_casadi_solver.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,14 @@ def test_model_solver(self):
model.events = [pybamm.Event("an event", var + 1)]
disc.process_model(model)
solver = pybamm.CasadiSolver(rtol=1e-8, atol=1e-8)
t_eval = np.linspace(0, 1, 100)
solution = solver.solve(model, t_eval)
np.testing.assert_array_equal(solution.t, t_eval)
np.testing.assert_array_almost_equal(
solution.y[0], np.exp(0.1 * solution.t), decimal=5
)

# Safe mode, without grid (enforce events that won't be triggered)
solver = pybamm.CasadiSolver(mode="safe without grid", rtol=1e-8, atol=1e-8)
solution = solver.solve(model, t_eval)
np.testing.assert_array_equal(solution.t, t_eval)
np.testing.assert_array_almost_equal(
Expand Down

0 comments on commit d24d4aa

Please sign in to comment.