refactor solver class

[martinjrobins]

see #16.

This class will take an instance of a model that has been processed by the Parameter (#18) and Spatial Discretisation (#20) classes. That is, the solver can assume that all Parameter, Variable and Operator nodes representing gradients (e.g. Laplace) nodes in the expression tree have been replaced with concrete Vector and Matrix nodes. That is, the expression tree now represents a concrete linear algebra expression

For now we only consider Models representing ODE equations and solvers implementing of methods of lines (MOL). But the implementation should be flexible enough to be extended to DAE equations, and time discretistion with a newton iteration solver rather than MOL.

The solver class needs to:

1. verify that the model rhs expressions are concrete linear algebra expressions
2. loop over time: 2.1. evaluate RHS as a function of current state variable y and t 2.2. step forward in time using RHS

Solver class should store the current state vector internally, and have an integrate function that can move the simulation forward in time, to allow for output/plotting over time, e.g.

model = ReactionDiffusionModel()
param_file = Parameters("file.dat")
param_file.process(model)
disc = FiniteVolumeDiscretisation(mesh_dx=0.1, something_else=2)
y0 = disc.process(model)
solver = RungeKuttaSolver(y0, tolerance=1e-5)
Tf = 1.0
nout = 100
nout_dt = Tf / nout
for i in range(nout):
   solver.integrate(nout_dt)
   plot(solver.y, disc)

[valentinsulzer]

Am having a look at this now. It's not obvious that the step-by-step integrate-and-plot is possible with adaptive time-steppers such as scipy.integrate.solve_ivp, so I will start by implementing something along the lines of

...
tspan = disc.mesh["time"]
solver = RungeKuttaSolver(tolerance = 1e-5)
solver.solve(model, y0, tspan) #y0 could be inside model, e.g. replace model.initial_conditions with y0

After solving, we want to have some objects that we can visualise (concentrations, electric potentials, current densities, fluxes, etc). Some of these will be explicit variables of the model (i.e. they have ODEs in model.rhs explicitly defined for them, e.g. concentration), while others will be derived variables (i.e. we can calculate them from explicit variables, e.g. flux). I'm not sure what solver.solve should set/return so that we can plot?

Suggestion:

1. define some variables for plotting in Model class:

   model.variables_for_plotting = {"c": c, "N": grad(c), etc} # values are pure expressions

2. set parameters and discretise model.variables_for_plotting as with other parts of the model
3. solve: solver sets a solution solver.t, solver.y
4. we can now evaluate model.variables_for_plotting and plot:

   solved_c = model.variables_for_plotting["c"].evaluate(solver.t, solver.y)
   plot(disc.mesh["whole_cell"].nodes, solved_c)

[martinjrobins]

I like the idea of replacing this model.initial_conditions with y0.

I don't think the model should hold a bunch of variables for plotting, as I'm sure you would get the case of two users wanting to use the same model, but also wanting entirely different plots. I think it would be pretty straightforward to plot an arbitrary derived list of derived variables, given a solution that consists of just y and t. you would just write something like:

to_plot = [c, grad(c), c**2, c*x]
to_plot = [disc.process_symbol(x).evaluate(t,y) for x in to_plot]
pybamm.plot(disc.mesh["whole_cell"].nodes, to_plot)

[valentinsulzer]

I agree that different users might want different plots, but I think the model should contain a dictionary {"variable name": variable} of all variables (explicit or derived), from which the user can choose which ones to plot.

For example, the user might just want to plot "the current" as defined by the model - but different models might have defined this differently (usually some functional of concentration and potentials).

Then, for example

to_plot = [model.variables[x] for x in ["concentration", "flux", "current"]] 
           + [model.variables["concentration"]**2]

model.variables can be discretised at the same stage as rhs and initial_conditions

[martinjrobins]

ok, yea, I can see that would be handy. So model.variables will be a dict of all the relevant variables (explicit and derived) for that particular model?