Casadi solver doesn't work with external variables

[TomTranter]

Summary
ODE solver allows you to specify external variables - it would be good if Casadi solver did to

At the very least there should be some error messages saying it's not implemented yet. If you try to run it you get runtime errors

RuntimeError: .../casadi/core/function_internal.hpp:1247: Input 1 (i1) has mismatching shape. Got 80-by-1. Allowed dimensions, in general, are:

• The input dimension N-by-M (here 81-by-1)
• A scalar, i.e. 1-by-1
• M-by-N if N=1 or M=1 (i.e. a transposed vector)
• N-by-M1 if K*M1=M for some K (argument repeated horizontally)
• N-by-P*M, indicating evaluation with multiple arguments (P must be a multiple of 1 for consistency with previous inputs)

This was produced with the following script

import pybamm
import numpy as np
model_options = {
    "thermal": "x-lumped",
    "external submodels": ["thermal"],
}
model = pybamm.lithium_ion.SPMe(model_options)
solver = pybamm.CasadiSolver()
sim = pybamm.Simulation(model, solver=solver)
t_eval = np.linspace(0, 0.01, 3)
T_av = 0.0
for i in np.arange(1, len(t_eval) - 1):
    dt = t_eval[i + 1] - t_eval[i]
    external_variables = {"X-averaged cell temperature": T_av}
    T_av += 1.0
    sim.step(dt, external_variables=external_variables)

[rtimms]

looks like this is because the padding of the statevector happens inside the SolverCallable class, but when using the Casadi solver the attributes self.casadi_rhs and self.casadi_algebraic use the casadi.Function objects directly, so the statevector doesn't get padded.

One option could be to pass an extra argument to the SolverCallable class telling it what kind of evaluate to use, e.g.

class Rhs(SolverCallable):
    "Returns information about rhs at time t and state y"

    def __init__(self, concatenated_rhs_fn, eval_type):
        self.concatenated_rhs_fn = concatenated_rhs_fn
        self.eval_type = eval_type

    def __call__(self, t, y):
        y = y[:, np.newaxis]
        y = add_external(y, self.y_pad, self.y_ext)
        if eval_type == "python":
            return self.concatenated_rhs_fn(t, y, self.inputs, known_evals={})[0][:, 0]
        elif eval_type == "casadi to python":
            return self.concatenated_rhs_fn(t, y, self.inputs_casadi).full()[:, 0]
        elif eval_type == "casadi":
            return self.concatenated_rhs_fn(t, y, self.inputs_casadi)   # not sure if this is exactly right for this one?

or maybe it is more efficient to create another class which pads then returns the casadi.Function directly (not sure if the if eval type is slow?).

[valentinsulzer]

I think the proper way of dealing with external variables / inputs in casadi is to pass them as "parameters" to the solver - otherwise there might be "variable is free" errors. This is why I didn't sort it out properly at the same time as the other solvers.

Separate but related - do you think external variables and inputs are fundamentally different objects which should be kept separate (as currently done) or should we unify the interface so that external variables are treated as an input that is not scalar?

[rtimms]

Ah ok - I wasn't too sure on how to correctly to pass to casadi, but it seems like you have a better understanding of this.

I don't think external variables and inputs are fundamentally different. Presumably we can keep all of the machinery for external variables in discretisation etc., but then pass the vectors in as inputs using the variable name? I think a unified interface would be cleaner.

[TomTranter]

They are pretty similar sounding. But I kind of understood it as inputs are changeable parameters that are never usually "solved-for" and external variables are fixed variables that are usually solved-for. This pretty much amounts to the same thing but I guess the analogy would be like linking model outputs to inputs in comsol. I'm not sure if we need to make the distinction clearer or just merge them. Would it simplify the underlying code?

[TomTranter]

Hey @tinosulzer is there a chance this can get done this week? Would really help me out. Cheers

[valentinsulzer]

ok will get the hacky solution done asap

[valentinsulzer]

Branch issue-775-casadi-external needs some work to get it fixed and merged into master, but you can at least merge it into your own branch for now and the test you posted above will work

[TomTranter]

@tinosulzer awesome thanks a lot. I merged it into my branch and it seems to be working. Do you have an idea for a better way of doing it that you are still working on. Want me to review when that's done or shall we get this signed off and merged in now?

[valentinsulzer]

Let's merge this in now, and I'll work on other issues as part of #784

[TomTranter]

Fine by me . I approved the PR