pybamm-team / pybamm

Fast and flexible physics-based battery models in Python

PyBaMM

PyBaMM (Python Battery Mathematical Modelling) is an open-source battery simulation package written in Python. Our mission is to accelerate battery modelling research by providing open-source tools for multi-institutional, interdisciplinary collaboration. Broadly, PyBaMM consists of (i) a framework for writing and solving systems of differential equations, (ii) a library of battery models and parameters, and (iii) specialized tools for simulating battery-specific experiments and visualizing the results. Together, these enable flexible model definitions and fast battery simulations, allowing users to explore the effect of different battery designs and modeling assumptions under a variety of operating scenarios.

PyBaMM uses an open governance model and is fiscally sponsored by NumFOCUS. Consider making a tax-deductible donation to help the project pay for developer time, professional services, travel, workshops, and a variety of other needs.

💻 Using PyBaMM

The easiest way to use PyBaMM is to run a 1C constant-current discharge with a model of your choice with all the default settings:

import pybamm

model = pybamm.lithium_ion.DFN()  # Doyle-Fuller-Newman model
sim = pybamm.Simulation(model)
sim.solve([0, 3600])  # solve for 1 hour
sim.plot()

or simulate an experiment such as a constant-current discharge followed by a constant-current-constant-voltage charge:

import pybamm

experiment = pybamm.Experiment(
    [
        (
            "Discharge at C/10 for 10 hours or until 3.3 V",
            "Rest for 1 hour",
            "Charge at 1 A until 4.1 V",
            "Hold at 4.1 V until 50 mA",
            "Rest for 1 hour",
        )
    ]
    * 3,
)
model = pybamm.lithium_ion.DFN()
sim = pybamm.Simulation(model, experiment=experiment, solver=pybamm.CasadiSolver())
sim.solve()
sim.plot()

However, much greater customisation is available. It is possible to change the physics, parameter values, geometry, submesh type, number of submesh points, methods for spatial discretisation and solver for integration (see DFN script or notebook).

For new users we recommend the Getting Started guides. These are intended to be very simple step-by-step guides to show the basic functionality of PyBaMM, and can either be downloaded and used locally, or used online through Google Colab.

Further details can be found in a number of detailed examples, hosted here on github. In addition, there is a full API documentation, hosted on Read The Docs. Additional supporting material can be found here.

Note that the examples on the default develop branch are tested on the latest develop commit. This may sometimes cause errors when running the examples on the pybamm pip package, which is synced to the main branch. You can switch to the main branch on github to see the version of the examples that is compatible with the latest pip release.

Versioning

PyBaMM makes releases every four months and we use CalVer, which means that the version number is YY.MM. The releases happen, approximately, at the end of January, May and September. There is no difference between releases that increment the year and releases that increment the month; in particular, releases that increment the month may introduce breaking changes. Breaking changes for each release are communicated via the CHANGELOG, and come with deprecation warnings or errors that are kept for at least one year (3 releases). If you find a breaking change that is not documented, or think it should be undone, please open an issue on GitHub.

🚀 Installing PyBaMM

PyBaMM is available on GNU/Linux, MacOS and Windows. We strongly recommend to install PyBaMM within a python virtual environment, in order not to alter any distribution python files. For instructions on how to create a virtual environment for PyBaMM, see the documentation.

Using pip

pip install pybamm

Using conda

PyBaMM is available as a conda package through the conda-forge channel.

conda install -c conda-forge pybamm

Optional solvers

The following solvers are optionally available:

• jax-based solver, see the documentation.

📖 Citing PyBaMM

If you use PyBaMM in your work, please cite our paper

Sulzer, V., Marquis, S. G., Timms, R., Robinson, M., & Chapman, S. J. (2021). Python Battery Mathematical Modelling (PyBaMM). Journal of Open Research Software, 9(1).

You can use the BibTeX

@article{Sulzer2021,
  title = {{Python Battery Mathematical Modelling (PyBaMM)}},
  author = {Sulzer, Valentin and Marquis, Scott G. and Timms, Robert and Robinson, Martin and Chapman, S. Jon},
  doi = {10.5334/jors.309},
  journal = {Journal of Open Research Software},
  publisher = {Software Sustainability Institute},
  volume = {9},
  number = {1},
  pages = {14},
  year = {2021}
}

We would be grateful if you could also cite the relevant papers. These will change depending on what models and solvers you use. To find out which papers you should cite, add the line

pybamm.print_citations()

to the end of your script. This will print BibTeX information to the terminal; passing a filename to print_citations will print the BibTeX information to the specified file instead. A list of all citations can also be found in the citations file. In particular, PyBaMM relies heavily on CasADi. See CONTRIBUTING.md for information on how to add your own citations when you contribute.

🛠️ Contributing to PyBaMM

If you'd like to help us develop PyBaMM by adding new methods, writing documentation, or fixing embarrassing bugs, please have a look at these guidelines first.

📫 Get in touch

For any questions, comments, suggestions or bug reports, please see the contact page.

📃 License

PyBaMM is fully open source. For more information about its license, see LICENSE.

✨ Contributors

This project follows the all-contributors specification. Contributions of any kind are welcome!

Click here to see a full list of our contributors' profiles.

pybamm's Issues

change boundary conditions format to {variable: {"left" => left_bc, "right" => right_bc}}

so, for example, you can set the boundary conditions in a model like so:

self.boundary_conditions = {N_e: {"left": 0, "right": 0}}

Add mesh and operators for spherical domain

Summary

• Add (dimensionless?) spherical domain to mesh.py, independently of the existing Cartesian domain
• Add grad and div operators (Finite Volumes) for the spherical domain

Implement Single Particle Model (SPM)

Summary
Implement the SPM. We require:

• a StandardParticle class which implements spherical diffusion on a spherically symmetric particle taking a boundary condition G for the flux on the particle surface.
• an SPM class which calls homogeneous_reactions to create G and then calls StandardParticle twice to create a negative electrode particle and a positive electrode particle.
• The SPM class should also contain an expression for the voltage in terms of the other fundamental variables (which are just c_n and c_p in this case). This will be contained within the variables dictionary of the model

convert expression tree to text

See #16.

It is generally useful to encode an expression tree to text, in order to:

1. identify common sub-trees
2. export to Fenics UML
3. export to Sympy expression
4. export to latex (note, once 3 is done a sympy expression can be converted to latex)
5. print out expression for the user (note, once 3 is done Sympy can pretty print an expression)

Add DAE option to solver

Add a new solver class DAESolver

• Import SUNDIALS (or other) in Dockerfile
• Tell the solver which equations are ODEs and which are algebraic; I think this should be done by the model class
• Write DAESolver class to solve the DAE system

add a dictionary attribute `variables` to model

After discussions in #21 , it would be useful for the BaseModel class to have an attribute variables with that maps "variable name" to a variable. This dict needs to be processed by parameters and discretisation just like rhs and initial_conditions. Will then be useful for visualisation.

Implement Doyle-Fuller-Newman (DFN) Model

Update
Scope of this issue has been extended to also implementing the SPM and SPMe by calling submodels since the forms of the models are so similar.

Summary
We need to implement a DFN class. To do this, we need:

• A butler_volmer function that creates a G given in terms of the other parameters and variables in the model
• An StefanMaxwellCurrent class which implements the conservation of current equations in the electrolyte as well as the McInnes equation for phi_e.
• A StandardElectrode class which implements current conservation and Ohm's law in the macroscale electrode
• A DFN class which:
  • makes a G using the butler_volmer function
  • creates the electrolyte model by calling StefanMaxwellDiffusion and StefanMaxwellCurrent
  • creates the macroscale electrode model by calling StandardElectrode.
  • create negative and positive electrode particles using StandardParticle. It must be somehow indicated that a particle must be inserted in every cell of the discretised model.
• Within DFN we should also add the following useful variables:
  • current
  • overpotentials
  • the terminal voltage

Docs build failing - guzzle_sphinx_theme

Doc build failing - see build report
Problem seems to be that guzzle_sphinx_theme is not being installed. It is in requirements-docs.txt, but not in set-up.py.

@martinjrobins is there a way to avoid having to remember to keep these two in sync (and is there anything else that also needs to be kept in sync?)

use an expression tree to encode a model

After discussions, we've decided to switch to using an expression tree to encode each model and submodel (see #14). The idea is that this would decouple the different aspects of PyBaMM (e.g. model/set_parmeters/discretisation/solver) and allow for a pipeline approach, where the expression tree is mutated by each aspect in turn.

This issue is to provide an initial implementation of the expression tree so that it can start being used by the rest of the code.

Possible class hierarchy:

1. class Symbol: base class, holds list of child nodes, has a (string) name
   • operators used to build up tree defined here (e.g. symbol + symbol returns Addition(symbol,symbol))
   • note: could we use anytree to provide basic tree functionality?: https://anytree.readthedocs.io/en/2.4.3/index.html
2. class UnaryOperator(Symbol): provides a base class for all the unary operators (e.g. Negate)
3. class BinaryOperator(Symbol): same but for binary ops
4. class Addition(BinaryOperator): implementation of the addition op,
5. class Laplacian(UnaryOperator): implementation of the laplacian op.
6. many other ops....
7. class Parameter(Symbol): a parameter (symbolic, no value). A Parameter has a list of domains (text) that it is valid over
8. class Variable(Symbol): a variable that will (eventually) be discretised over the mesh. At the moment it is purely symbolic. A variable has a list of domains (text) that it is valid over
9. class Vector(Symbol): a node representing a vector (the spatial discretisation class/function will add these in)
10. class Matrix(Symbol): a node representing a matrix (the spatial discretisation class/function will add these in)
11. class Scalar(Symbol): a node representing a scalar (the set_parameters function will add these in place of Parameter). A Scalar has a list of index ranges, and a value for each range
12. class IndependentVariable(Symbol): a node representing an independent variable (space or time)

iteration through tree

provide functionality for iterating through the tree, e.g.

[node.name for node in PreOrderIter(f)] # pre-order traversal
[node.name for node in PostOrderIter(f)] # post-order traversal
[node.name for node in LevelOrderIter(f)] # level-order (i.e. breadth-first)

(note: these examples come from the anytree library)

other

there will be other functionality that we will need. e.g converting to/from other types of expression trees like Fenics UML, but this issue is for the initial implementation, other things can go in separate issues...

Representation of expression tree classes

Currently the only class with __repr__ implemented is Symbol. Other classes need a more accurate __repr__ than just inheriting from Symbol.

Ensure h**2 convergence even with non-uniform mesh

Finite Volumes convergence is not guaranteed to be h**2 if the grid is not uniform
(Diskin, Boris, and James L. Thomas. "Notes on accuracy of finite-volume discretization schemes on irregular grids." Applied numerical mathematics 60.3 (2010): 224-226.)

To do: Change Mesh class so that convergence is h**2 even if the mesh is not uniform

IndependentVariable: a node representing an independent variable (space or time)

Add to expression tree hierarchy: class IndependentVariable(Symbol): a node representing an independent variable (space or time)

When discretising, replace model.initial_conditions with concatenated `y0`

Small fix to issue #20 after discussions in #21 - discretisation should concatenate model.initial_conditions in place rather than returning y0.

Increase Modularity of Models

I have pushed a branch: component_class to my fork of PyBaMM if you want an example. On that branch, I have briefly explore this potential class structure for models.

The basic idea is to have two sets of classes: Models and SubModels each of which have common methods and variables. Then we derive other classes from these such as:
class SPM(Model):
which then inherits the methods common to Model. It seems that this would make it easier to add a new type of model as we inherit all the standard methods.

The Models are a two layer structure with

• Model
  • SPM
  • SPMe,
  • DFN,
  • LeadAcid,
  • User
  • etc

The superclass 'Model' defines the methods and variables common to all models such as .initial_conditions() and .rhs()
The subclasses on level 2. then just define things specific to that model, e.g the domain and the type of submodel to use on that domain.

SubModels are arranged in a three layer structure:

• SubModel
  • Particle
    • StandardParticle
    • DeformingParticle
    • PhaseFieldParticle
  • Electrolyte
    • NernstPlank,
    • StefanMaxwell,
  • Electrode
    • StandardElectrode
    • BinderElectrode
  • CurrentCollector
    • StandardCurrentCollector

All SubModels should have a common structure so that they can be interacted with by Models in a consistent way. Therefore a superclass SubModel that contains common methods and variables is used (might just be conceptually useful as, for example, .rhs() will be different for each submodel).

We then have a different type of SubModel for each sub_domain: Particle, Electrolyte, etc. These define the methods and variables common to each of their domains such as setting initial conditions but allow for further specification. Finally, on level 3. a particular SubModel is choosen for a particular domain. Here specifics like the exact model equations on that domain are stated.

Motivation
Increase modularity of the model section of the code

add a `Function` node to expression tree

for Parameters that vary with time, user puts in the name of the function (functions defined in .py in the same directory as the parameter file). Then the parameter class replaces Parameter with Function of that name

discritisation of boundary conditions

I was just looking at the discretisation of the boundary conditions. At the moment this seems to be mainly done in the gradient and divergence functions of MatrixVectorDiscretisation, which both do this:

       # Add boundary conditions if defined
        if symbol.id in boundary_conditions:
            lbc, rbc = boundary_conditions[symbol.id]
            discretised_symbol = self.concatenate(lbc, discretised_symbol, rbc)
        gradient_matrix = self.gradient_matrix(domain)
        return gradient_matrix * discretised_symbol

So it looks like extra ghost nodes are added to the current state vector to implement the bc. Should this be done instead when the pybamm.Variable node is processed in process_symbol? Also, the way the boundary conditions are implemented would depend on the method, so this should probably be done in FiniteVolumeDiscretisation? @tinosulzer?

Mesh class

Create a mesh class that takes a Parameter class with mesh parameters, for eventual input to the spatial discretisation class (#20)

eg.

model = ReactionDiffusionModel()
param_file = Parameters("file.dat")
param_file.process(model)
mesh = Mesh(param_file)
disc = FiniteVolumeDiscretisation(mesh, something_else=2)
y0 = disc.process(model)

Update sphinx setup for online documentation

• Configure sphinx to update with new documentation
• Write how this was done for future reference

refactor Simulation class

Refactor Simulation class to run the whole process (create model, set parameters, discretise, solve, visualise). Needs to take optional inputs for the parameter class, parameter file, mesh/discretisation, solver, and set defaults if these are None. Can later add capability to save solution, load a saved solution, etc.

class Simulation(object):
   def __init__(
        self,
        model,
        parameter_values=None,
        discretisation=None,
        solver=None,
        name="unnamed",
    ):
        # Defaults
        if parameter_values is None:
           # default parameter_values
        if discretisation is None:
           # default discretisation
        if solver is None:
            # default solver

        self.model = model
        self.parameter_values = parameter_values
        self.discretisation = discretisation
        self.solver = solver
        self.name = name

    def set_parameters(self):
        self.parameter_values.process(self.model)

    def discretise(self):
        y0 = self.discretisation.process(self.model)
        return y0

    def solve(self, y0):
        self.solver.solve(self.model, y0)

    def run(self):
        self.set_parameters()
        y0 = self.discretise()
        self.solve(y0)

    def load(self):

    def save(self):

    def plot(self):

documentation formatting for expression tree nodes

The formatting for the documentation on the expression tree node classes does not match the rest of pybamm, need to make these consistent

Explore FEniCS integration

Summary
Add functionality that allows the existing models to be solved using FEniCS instead of Finite Volumes

Motivation
Ideally we should be able to switch between different spatial discretisations (e.g. Finite Volumes and Finite Elements) without affecting the rest of the code (especially the model)

Additional context
Non-exhaustive to-do list:

• Add a FEniCS mesh to Mesh()
• Add an option to return a weak form?
• Implement a time integrator that allows us to use FEniCS

Process scalar initial conditions

Process initial conditions to take in combinations of scalars and return a vector with the right shape and value.

• Check that the initial condition is a combination of scalars (evaluates to a numbers.Number), if not raise NotImplementedError
• Change the function scalar_to_vector to vector_of_ones that takes in a domain and returns a Vector of the right shape with value 1.

visualise simulation results

Add class(es)/function(s) to visualise the result of a simulation. Options:

• time vs variable plot for variables that are only functions of time (current, voltage, temperature, etc)
• Slider plot (space vs variable at a time definable by a slider) for variables that depend on space and time (concentration, etc)
• Export for plotting using external software (e.g. pgfplots, Paraview)
• others?

Domains as strings

Describe the bug
The Domain class checks that domain is an iterable and raises a TypeError if it isn't.
However, a str domain is not caught by this as iter("a string") == True.
This isn't caught by the test in TestIndependentVariable; test should be

with self.assertRaises(TypeError):
   pybamm.IndependentVariable("a", domain="test")

Suggested fix
If domain is a string, replace with list(domain).

components.py

• is there any way to avoid passing grad and div into these functions. the reason being that if we are doing particles and y-z for the current collectors then we
  would have to pass in grad_x, div_x, grad_y, div_y, grad_r, ... etc into model_class. Obviously we could group them together but if we could just c.div_x for variables
  in the x-domain it may be cleaner?
• should we to pass a vector of ones and zeros to indicate ODEs or DAEs at this point

Add lithium ion parameters

Summary
Edit parameters.py to load lithium-ion instead of lead-acid

Motivation
We should be able to switch between lithium-ion and lead-acid parameters simply by changing the name of the input file

To do

• Create defaults_lithium-ion.csv and add lithium-ion parameters
• Edit parameters.py as appropriate - e.g. with any new dimensionless parameters

overloaded binary operators on `Symbol` take numbers

see #36. It might be useful for the overloaded binary operators onSymbol (e.g. __mul__) be able to take Numbers (e.g. int'). They would then implicitly convert this to a Scalar` node.

Note: also need to add __pow__ to this list.

refactor spatial discretisation class

see #16.

the new spatial discretisation base class and derived classes need to process a Model, along with an input mesh in order to:

1. count the number of variables needed, and therefore the size of the state vector
2. replace all Variable nodes with Vector nodes
3. replace all spatial gradient operators with Matrix nodes

example use:

model = ReactionDiffusionModel()
param_file = Parameters("file.dat")
param_file.process(model)
disc = FiniteVolumeDiscretisation(mesh_dx=0.1, something_else=2)
y0 = disc.process(model)

My current thoughts are that the discretisation class should also generate initial conditions (perhaps from functions contained in the Variable nodes of the expression tree?), and output the initial state vector from the process function so that it can be used by the solver class

only do discretisation of gradient and divergence once

for an expression with more than one grad(c) or div(c), the BaseDiscretisation will discretise each of them individually. We should identify these sub-expressions and replace them all with the same Matrix

note that identifying sub-expressions is also in #55, so this will overlap a bit

testing TODO

• general convergence test
• test convergence of complicated model to simple
• analytical solutions to specific models
• method of manufactured solutions

Implement the Single Particle Model with Electrolyte (SPMe)

Note: #72 Should be completed before this issue.

Summary
We require an SPMe class. This class builds upon the SPM class by adding an electrolyte model and modifying the voltage expression. This will be done by calling the StefanMaxwellDiffusion class to add a model for diffusion in the electrolyte. The StefanMaxwellDiffusion class will require a G as an input. You can either call the homogenous_reactions function again to create this G or it can be somehow extracted from SPM. The voltage expression should just be overwritten.

Note: SPMe could also be implemented without calling SPM and instead build the same things inside SPMe again but this might be a good exercise to see how easy it is to add physics to a model.

Refactor HomogenousReactions class

This is a pre-requisite for #17

Refactor this class to that it simply returns an symbolic expression for the Homogeneous reactions.

Should this really be a class? A function is probably sufficient

refactor Parameters to process an expression tree

See #16.

After a model has been created (see #19 and #17), the set_parameters() function takes an expression tree, and replaces all instances of Parameter nodes with Scalar nodes. e.g.

model = ReactionDiffusionModel()
param_file = Parameters("file.dat")
param_file.process(model)

add a list of index ranges to `Scalar`, and a value for each range

"A Scalar has a list of index ranges, and a value for each range". This is not implemented

optimise adjacent matricies and vectors after discretisation

After the discretisation step, there is the possibility to optimise adjacent Matrix, Vector and Scalar nodes, like in #54.

This might be best to do just once after the discretisation step, which would also solve #54.

Combining submodels

Following from #49

If each submodel (e.g. StefanMaxwellDiffusion) is a Model in its own right, we should have a function which allows for the simple concatenation of two submodels e.g. the functionality:

basemodel.create_from_submodels([submodelA, submodelB, submodelC, etc])

This would concatenate the rhs, initial_conditions, boundary_conditions, variables dictionaries.

It is possible that overlapping variable names may become an issue (e.g. concentration in negative and positive particle if within the particle expression they are simply called c). By associating a unique name to each submodel when it is initially defined, this could be avoided.

Process parameters for boundary conditions

In ParameterValues.process_model, the processing of boundary conditions assumes that boundary conditions are a single equation, not a dictionary ({"left": , "right": }).

Also, ParameterValues.process_model is currently untested.

ParameterValues class for a specific model

Building on issue #18 / PR #35

For a particular, model, parameter values may be needed that are combinations of "raw" parameter values. For example, the total electrode length L = Ln + Ls + Lp, or any dimensionless parameter.

It doesn't make sense to write these all in the csv file that gets imported, so there needs to be a class that calculates the "derived" parameters, maybe split by domain ("negative_electrode", "separator", "positive_electrode" or "whole_cell"). There may be overlap between the domains.
This class should also contain scales for non/re-dimensionalisation

For example,

class DiffusionParameterValues(BaseParameters):
    def __init__(self, base_parameters={}, optional_parameters={}):
       super().__init__(base_parameters, optional_parameters)
       
       self.scales = {"length": self.raw["Ln"] + self.raw["Ls"] + self.raw["Lp"],
                      "time": ...,
                      etc}
       self.whole_cell = {"L": self.raw["Ln"] + self.raw["Ls"] + self.raw["Lp"],
                          "Cd": self.scales["length"]**2/self.raw["D"]/self.scales["time"],
                          etc}
       self.negative_electrode = {"L": self.raw["Ln"], etc}
       self.separator = {"L": self.raw["Ls"], etc}
       self.positive_electrode = {"L": self.raw["Ln"], etc}

Maybe make DiffusionParameterValues emulate a dictionary, similarly to #19 , with keys being domain names.

Now if we want to include new parameters, say for a temperature model, we will need more scales and derived parameters, but these should be in a separate class as they might call upon raw values that would not necessarily be given to a DiffusionParameterValues class.

class TemperatureParameterValues(BaseParameters):
    def __init__(self, base_parameters={}, optional_parameters={}):
       super().__init__(base_parameters, optional_parameters)
       
       self.scales = {"Temp": self.raw["T_init"]}
       self.whole_cell = {"k": self.raw["D_temp"]/ bla bla}
       self.negative_electrode = {}
       self.separator = {}
       self.positive_electrode = {}

Then a BatteryParameterValues could multiply-inherit from DiffusionParameterValues and TemperatureParameterValues, combining the dictionary attributes from each of these (not sure how to do this)

@martinjrobins @Scottmar93 do you think this is a good approach? See previous implementation (quite different) for an idea of which derived parameters need to be implemented.

Non-spatial unary operators

Symbol types might need __neg__, __pow__, __abs__, others?

Using Binary Operations

Getting issues when trying to write

c = Variable('c') 
 1 - c  

and also

 -c 

Does the BinaryOperator class account for the case of integer inputs or no input on one side but symbols on the other?

Bug in parameter_values.process_symbols

Describe the bug
Binary operators are not being processed properly by ParameterValues. At line 119,

new_left = self.process_symbol(symbol.children[0])
new_right = self.process_symbol(symbol.children[1])

the children are being switched by the first line, so new_left is processed again by the second line:

# symbol.children = (left, right)
new_left = self.process_symbol(symbol.children[0])
# symbol.children = (right, new_left)
new_right = self.process_symbol(symbol.children[1])
# symbol.children = (right, new_new_left)

Solution
Possibly a hacky solution, but we can get around this by pre-setting pointers to left and right

left, right = symbol.children
new_left = self.process_symbol(left)
# symbol.children = (right, new_left)
new_right = self.process_symbol(right)
# symbol.children = (new_left, new_right)

refactor solver class

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)

optimise adjacent scalars in tree after set parameters

After ParameterValues swaps out all the Parameter for Scalar nodes, there is the possibility that you get things like Scalar(4) * Scalar(5) in the tree, which can be combined into Scalar(20)

refactor BaseModel to hold an expression tree

See #16.

The base model will internally hold an expression tree representing the rhs of an ode model. e.g.

class BaseModel:
   def __init__(self):
         self._rhs = ... # lhs is a dictionary mapping variables to expressions
         self._initial_conditions = ... # a dictionary mapping expression to expressions representing the initial conditions
         self._boundary_conditions = ... # a dictionary mapping expression to expressions representing the boundary conditions
   @property
   def rhs(self):
          return self._rhs
   ...

Note: rather than having rhs return a dict, it might be nicer to have BaseModel emulate a dictionary, so you could write model[variable] and get the rhs expression associated with that variable (see https://docs.python.org/3/reference/datamodel.html?emulating-container-types#emulating-container-types)

Add a solver test to StefanMaxwellDiffusion tests

Summary

Add a solver test to StefanMaxwellDiffusion tests.

Check convergence of Method of Lines

Using method of lines:

• Finite Volumes spatial discretisation -> ODEs
• scipy.integrate.solve_ivp(..., method='BDF') to solve ODEs

How to check convergence to a manufactured solution, both in space and time?
Would want to do this in test_simulation.py

Refactor StefanMaxwellDiffusion class

This issue must be completed prior to refactoring the reaction diffusion model #17

StefanMaxwellDiffusion should inherit from BaseModel so that it is a complete model in its own right. It will contain a dictionary with the set of equations in the electrolyte implemented in symbolic form employing the expression tree.

Will looks something like:

class StefanMaxwellDiffusion(BaseModel): 
        def __init__(self): 
            self._rhs = {c: "epression"}
            self._initial_conditions = .... etc 

           self._variable = {"c": c, ...} 

EDIT.
The keys should be of class Symbol expect for the _variable dictionary.

discretisation of different parameter values in different domains

Different parameters can have different values in different domains. This can be represented in a model by a concatentation of Parameters, which will change to a concat of Scalars. The discretisation class need to change this to a concatenation of Vectors.

class StateVector for expression tree

Need a node representing the state vector, to replace Variable in the discretisation step.

Should have an internal slice object (https://docs.python.org/dev/library/functions.html#slice) that can be None if the node returns the whole state vector

for example, from @tinosulzer:

class VariableVector(pybamm.Symbol):
    """A vector that depends on an input y."""

    def __init__(self, y_slice, name=None, parent=None):
        if name is None:
            name = str(y_slice)
        super().__init__(name=name, parent=parent)
        self._y_slice = y_slice

    @property
    def y_slice(self):
        return self._y_slice

    def evaluate(self, y):
        return y[self._y_slice]

refactor reaction diffusion model

See #16.

refactor the ReactionDiffusionModel class to inherit from BaseModel (see #19), implementing a reaction diffusion model