Look into using smooth approximations for min, max, heaviside

[valentinsulzer]

Functions with discontinuous derivatives can cause problems for adaptive time-stepping algorithms: https://discourse.julialang.org/t/handling-instability-when-solving-ode-problems/9019/5

Need to look into whether this is actually a problem for the solvers in pybamm. I suspect it might be given that linear interpolants have been giving similar issues.

If so, we could add a setting where pybamm automatically provides the smooth form instead of the exact form. e.g. by setting pybamm.settings.smooth_approximation_steepness, which is k in the link above, and setting k to zero could give the exact form.

[valentinsulzer]

Using smooth approximations definitely helps in a toy model:

import pybamm
def smooth_heaviside(left, right, k):
    return (1 + pybamm.tanh(k * (right - left))) / 2


def rhs_exact(x):
    return (x <= 1) * x + (x > 1) * (2 * x - 1)


def rhs_smooth(x, k):
    return smooth_heaviside(x, 1, k) * x + smooth_heaviside(1, x, k) * (2 * x - 1)
    # return (v2 - v1) * pybamm.tanh(k * (x - switch)) / 2 + (v1 + v2) / 2


x = pybamm.Variable("x")
y = pybamm.Variable("y")

model_exact = pybamm.BaseModel(name="exact")
model_exact.rhs = {x: rhs_exact(x)}
model_exact.initial_conditions = {x: 0.5}
model_exact.variables = {"x": x, "rhs": rhs_exact(x)}

model_smooth = pybamm.BaseModel(name="smooth")
k = pybamm.InputParameter("k")
model_smooth.rhs = {x: rhs_smooth(x, k)}
model_smooth.initial_conditions = {x: 0.5}
model_smooth.variables = {"x": x, "rhs": rhs_smooth(x, k)}

pybamm.set_logging_level("INFO")
# solver = pybamm.ScikitsDaeSolver()
solver = pybamm.CasadiSolver(extra_options_setup={"print_stats": True})
print("Exact-------------------------")
sols = [solver.solve(model_exact, [0, 2])]
for k in [1, 10, 50]:
    print(f"Smooth, k={k}-------------------------")
    # solver = pybamm.ScikitsDaeSolver()
    solver = pybamm.CasadiSolver(extra_options_setup={"print_stats": True})
    sol = solver.solve(model_smooth, [0, 2], inputs={"k": k})
    sols.append(sol)

pybamm.dynamic_plot(sols, ["x", "rhs"])

Results:

Exact: 4ms
Smooth (independent of k): 1ms
Solutions line up well for k=10,100 (k=1 is too small)

Not super obvious from the stats printed by casadi where the improvement comes from though