pybamm-team/liionpack

A battery pack simulation tool that uses the PyBaMM framework

liionpack's Issues

Add contributing and code of conduct

We can take the PyBaMM ones as an example.

pack voltage bug

Event handling

Currently there is none except monitoring the voltage limits on all cells manually and stopping the simulation. Would be amazing if we can leverage some of the event handling in PyBaMM

Acknowledge funding

Possible optimization in `solve_circuit`

If I understand correctly, solve_circuit constructs the A and z matrices each time it is called, then solves Ax=z. Doing some very basic profiling, around 85% of the time is spent setting up the A matrix and 15% solving the linear system (this is for the load_netlist.py example, I haven't tested how this scales with pack size).

Does the A matrix change between iterations or is it always the same? If it always stays the same, constructing it once outside of solve_circuit would significantly speed up solve_circuit, at least for the circuit size I tried

Use Dask to run the cell simulations in parallel

This isn't really an issue but more of a feature request or enhancement. Would it be possible to use Dask to run the cell simulations in parallel? I created a basic example (see below) of running the SPMe model for several discharges in parallel using Dask. Compare elapsed time with and without Dask by commenting out the appropriate section in main(). Elapsed times are given in the table below when running on an 8-core CPU.

Dask is made for massive parallelization and it's fairly easy to setup for CPU and GPU clusters. If it can be used with liionpack then it could provide a huge performance boost for large pack simulations. I haven't tried Casadi's parallel features beyond running on a single CPU but I don't think it will be easy to scale compared to using Dask.

Example	Elapsed time
No Dask	8.57 seconds
Dask	3.83 seconds

import matplotlib.pyplot as plt
import pybamm
import time
from dask.distributed import Client

def generate_plots(discharge, t, capacity, current, voltage):

    def styleplot(ax):
        ax.legend(loc='best')
        ax.grid(color='0.9')
        ax.set_frame_on(False)
        ax.tick_params(color='0.9')

    _, ax = plt.subplots(tight_layout=True)
    for i in range(len(discharge)):
        ax.plot(t[i], current[i], label=f'{discharge[i]} A')
    ax.set_xlabel('Time [s]')
    ax.set_ylabel('Current [A]')
    styleplot(ax)

    _, ax = plt.subplots(tight_layout=True)
    for i in range(len(discharge)):
        ax.plot(t[i], voltage[i], label=f'{discharge[i]} A')
    ax.set_xlabel('Time [s]')
    ax.set_ylabel('Terminal voltage [V]')
    styleplot(ax)

    _, ax = plt.subplots(tight_layout=True)
    for i in range(len(discharge)):
        ax.plot(capacity[i], voltage[i], label=f'{discharge[i]} A')
    ax.set_xlabel('Discharge capacity [Ah]')
    ax.set_ylabel('Terminal voltage [V]')
    styleplot(ax)

    plt.show()

def run_simulation(dis, t_eval):

    model = pybamm.lithium_ion.SPMe()

    param = model.default_parameter_values
    param['Current function [A]'] = '[input]'

    sim = pybamm.Simulation(model, parameter_values=param)
    sim.solve(t_eval, inputs={'Current function [A]': dis})

    return sim.solution

def main(client):
    tic = time.perf_counter()

    discharge = [4, 3.5, 3, 2.5, 2, 1.8, 1.5, 1]  # discharge currents [A]
    t_eval = [0, 4000]                            # evaluation time [s]

    # No Dask
    # ------------------------------------------------------------------------

    # label = 'no Dask'

    # sols = []
    # for dis in discharge:
    #     sol = run_simulation(dis, t_eval)
    #     sols.append(sol)

    # Dask
    # ------------------------------------------------------------------------

    label = 'Dask'

    lazy_sols = client.map(run_simulation, discharge, t_eval=t_eval)
    sols = client.gather(lazy_sols)

    # ------------------------------------------------------------------------

    t = []
    capacity = []
    current = []
    voltage = []

    for sol in sols:
        t.append(sol['Time [s]'].entries)
        capacity.append(sol['Discharge capacity [A.h]'].entries)
        current.append(sol['Current [A]'].entries)
        voltage.append(sol["Terminal voltage [V]"].entries)

    toc = time.perf_counter()
    print(f'Elapsed time ({label}) = {toc - tic:.2f} s')

    generate_plots(discharge, t, capacity, current, voltage)

if __name__ == '__main__':
    client = Client()
    print(client)
    main(client)

Array split does not result in equal division

Splitting the arrays for attributes self.split_index and self.htc in solvers.py throws an error if the number of cells is not evenly divisible by the number of processes nproc. For example, using Np = 32, Ns = 20 with nproc = 6 does not work for the Dask and Ray solvers. This can be avoided using np.array_split but this creates another error about a singular matrix caused by divide by zero in netlist_utils.py for g[R_map] = 1 / value[R_map].

Switch to object-oriented programming?

This would be slightly more "pythonic". e.g.

• A Circuit class to handle loading a netlist, drawing, and solving
• A Solver class
• A PackSolution class to store, plot, save the solution

Though of course, there is nothing wrong with function-oriented programming as the code is currently designed.

Add packages to install requires

The following packages are not installed as requirements:

• plotly
• tqdm

Use lipack-specific logger

To enable setting a different logging level from pybamm.
We can just copy over the logger.py code from pybamm

Add coordinates to generated netlists

For scatter plot the coordinates of the cells are needed. These can be automatically generated for built netlists.

Change file defaults

Some functions use defaults for filenames and directories which are specific for a local installation (e.g. load_data.py l11, or batteryrig_utils.py l164). We should either remove the defaults or include the default dataset in the repo as an example.

Better initial guess for current

This is currently set by the netlist which can be created with any arbitrary current. It should come from the first step of the experiment

draw_circuit no longer showing values

Make initial_soc work here instead of using init_funcs

Make latex figures work in tests

In Github actions

======================================================================
ERROR: test_draw_circuit (tests.unit.test_plots.plotsTest)

Traceback (most recent call last):
File "/home/runner/work/liionpack/liionpack/tests/unit/test_plots.py", line 14, in test_draw_circuit
lp.draw_circuit(net)
File "/home/runner/work/liionpack/liionpack/liionpack/plots.py", line 54, in draw_circuit
cct.draw(**kwargs)
File "/opt/hostedtoolcache/Python/3.9.7/x64/lib/python3.9/site-packages/lcapy/netlistmixin.py", line 1636, in draw
return cct.sch.draw(filename=filename, **kwargs)
File "/opt/hostedtoolcache/Python/3.9.7/x64/lib/python3.9/site-packages/lcapy/schematic.py", line 784, in draw
self.tikz_draw(filename=filename,
File "/opt/hostedtoolcache/Python/3.9.7/x64/lib/python3.9/site-packages/lcapy/schematic.py", line 654, in tikz_draw
pdfconverter.to_png(pdf_filename, root + '.png', self.dpi)
File "/opt/hostedtoolcache/Python/3.9.7/x64/lib/python3.9/site-packages/lcapy/system.py", line 163, in to_png
raise RuntimeError("""Could not convert pdf to png, tried: %s. Check that one of these programs is installed.""" % ', '.join([conversion for conversion, func in conversions]))
RuntimeError: Could not convert pdf to png, tried: ghostscript, convert, pdftoppm. Check that one of these programs is installed.

Improve output data format and include cell currents by default

Unused functions in utils module

The following functions are not being used by the simulation and some only appear in the unit tests:

• interp_current() defined in liionpack/utils.py and used in tests/unit/test_utils.py
• _z_from_plane() defined in liionpack/utils.py
• _fit_plane() defined in liionpack/utils.py
• read_cfd_data() defined in liionpack/utils.py and used in tests/unit/test_utils.py
• get_linear_htc() defined in liionpack/utils.py and used in tests/unit/test_utils.py
• get_interpolated_htc() defined in liionpack/utils.py and used in tests/unit/test_utils.py

Is there a reason why these functions are defined in liionpack/utils.py? Can they be removed from the package?

Investigate PySpice support

HELP readthedocs nightmare!

I have no idea why but it seems like readthedocs is stuck reading from an old commit. It won't update the additional pages. Seems to update index ok. I'm using MKDocs which is supposed to make life simple. If I start a server locally it all looks ok :(

Casadi solver does not scale with CPU cores

I ran a 32 battery cell example (Np = 16, Ns = 2) on a 64 CPU core workstation where I adjusted the nproc parameter from 2 to 64. Based on the results shown below, it appears that the Casadi solver does not scale beyond 16 CPU cores. I get similar elapsed times for anything above 16 cores. I'm not sure how the Casadi map function executes in parallel so I have no idea how to optimize the problem to take advantage of many CPUs. Has anyone used liionpack with more than 16 CPU cores and did you notice any speed improvements above 16 cores?

Fix benchmark workflows

Fix the errors in the existing benchmark workflows

saving output for a simulation reaching high/low voltage cut-off

If a cell in the simulation reaches high/low cut-off voltage, the output is not saved. The error looks something like the following:
**Solving Pack: 25%|████████████████ | 156/630 [00:04<00:13, 35.83it/s]High voltage limit reached
Solving Pack: 25%|████████████████▍ | 159/630 [00:04<00:13, 34.70it/s]
Traceback (most recent call last):
File "/Users/uvk/Desktop/BatterySimulations/TestLiionpsck/liionpackNov8/BatSim.py", line 39, in
lp.plot_output(output)
File "/Users/uvk/Desktop/BatterySimulations/TestLiionpsck/liionpackNov8/liionpack/plots.py", line 282, in plot_output
plot_pack(output)
File "/Users/uvk/Desktop/BatterySimulations/TestLiionpsck/liionpackNov8/liionpack/plots.py", line 237, in plot_pack
ax2.plot(time, i_pack, color="blue", label="simulation")
File "/Users/uvk/anaconda2/envs/py3/lib/python3.9/site-packages/matplotlib/axes/_axes.py", line 1605, in plot
lines = [*self._get_lines(*args, data=data, kwargs)]
File "/Users/uvk/anaconda2/envs/py3/lib/python3.9/site-packages/matplotlib/axes/_base.py", line 315, in call
yield from self._plot_args(this, kwargs)
File "/Users/uvk/anaconda2/envs/py3/lib/python3.9/site-packages/matplotlib/axes/_base.py", line 501, in _plot_args
raise ValueError(f"x and y must have same first dimension, but "
ValueError: x and y must have same first dimension, but have shapes (159,) and (160,)

So i changed the solver_utils.py, where the time,current,voltage needs to be added to list before checking for the limits of the high/low cut-off voltages. Now I can see the output of the simulation. Here is the modified code:

#
# Solver utilities
#

import casadi
import pybamm
import numpy as np
import time as ticker
import liionpack as lp
from tqdm import tqdm


def _mapped_step(model, solutions, inputs_dict, integrator, variables, t_eval):
    """
    Internal function to process the model for one timestep in a mapped way.
    Mapped versions of the integrator and variables functions should already
    have been made.
    Parameters
    ----------
    model : :class:`pybamm.lithium_ion.BaseModel`
        The built battery model
    solutions : list of :class:`pybamm.Solution` objects for each battery
        Used to get the last state of the system and use as x0 and z0 for the
        casadi integrator
    inputs_dict : list of inputs_dict objects for each battery
    integrator : mapped casadi.integrator
        Produced by `_create_casadi_objects`
    variables : mapped variables evaluator
        Produced by `_create_casadi_objects`
    t_eval : float array of times to evaluate
        Produced by `_create_casadi_objects`
    Returns
    -------
    sol : list
        Solutions that have been stepped forward by one timestep
    var_eval : list
        Evaluated variables for final state of system
    """
    len_rhs = model.concatenated_rhs.size
    N = len(solutions)
    if solutions[0] is None:
        # First pass
        x0 = casadi.horzcat(*[model.y0[:len_rhs] for i in range(N)])
        z0 = casadi.horzcat(*[model.y0[len_rhs:] for i in range(N)])
    else:
        x0 = casadi.horzcat(*[sol.y[:len_rhs, -1] for sol in solutions])
        z0 = casadi.horzcat(*[sol.y[len_rhs:, -1] for sol in solutions])
    # t_min = [0.0]*N
    t_min = 0.0
    inputs = []
    for temp in inputs_dict:
        inputs.append(casadi.vertcat(*[x for x in temp.values()] + [t_min]))
    ninputs = len(temp.values())
    inputs = casadi.horzcat(*inputs)
    # p = casadi.horzcat(*zip(inputs, external_variables, [t_min]*N))
    # inputs_with_tmin = casadi.vertcat(inputs, np.asarray(t_min))
    # Call the integrator once, with the grid
    timer = pybamm.Timer()
    tic = timer.time()
    casadi_sol = integrator(x0=x0, z0=z0, p=inputs)
    integration_time = timer.time()
    nt = len(t_eval)
    xf = casadi_sol["xf"]
    # zf = casadi_sol["zf"]
    sol = []
    xend = []
    for i in range(N):
        start = i * nt
        y_sol = xf[:, start:start + nt]
        xend.append(y_sol[:, -1])
        # Not sure how to index into zf - need an example
        sol.append(pybamm.Solution(t_eval, y_sol, model, inputs_dict[i]))
        sol[-1].integration_time = integration_time
    toc = timer.time()
    lp.logger.debug(f"Mapped step completed in {toc - tic}")
    xend = casadi.horzcat(*xend)
    var_eval = variables(0, xend, 0, inputs[0:ninputs, :])
    return sol, var_eval


def _create_casadi_objects(I_init, htc, sim, dt, Nspm, nproc, variable_names):
    """
    Internal function to produce the casadi objects in their mapped form for
    parallel evaluation
    Parameters
    ----------
    I_init : float
        initial guess for current of a battery (not used for simulation).
    htc : float
        initial guess for htc of a battery (not used for simulation).
    sim : :class:`pybamm.Simulation`
        A PyBaMM simulation object that contains the model, parameter values,
        solver, solution etc.
    dt : float
        The time interval (in seconds) for a single timestep. Fixed throughout
        the simulation
    Nspm : int
        Number of individual batteries in the pack.
    nproc : int
        Number of parallel processes to map to.
    variable_names : list
        Variables to evaluate during solve. Must be a valid key in the
        model.variables
    Returns
    -------
    integrator : mapped casadi.integrator
        Solves an initial value problem (IVP) coupled to a terminal value
        problem with differential equation given as an implicit ODE coupled
        to an algebraic equation and a set of quadratures
    variables_fn : mapped variables evaluator
        evaluates the simulation and output variables. see casadi function
    t_eval : float array of times to evaluate
        times to evaluate in a single step, starting at zero for each step
    """
    inputs = {
        "Current function [A]": I_init,
        "Total heat transfer coefficient [W.m-2.K-1]": htc,
    }
    solver = sim.solver

    # Initial solution - this builds the model behind the scenes
    # solve model for 1 second to initialise the circuit
    t_eval = np.linspace(0, 1, 2)
    sim.solve(t_eval, inputs=inputs)

    # Step model forward dt seconds
    t_eval = np.linspace(0, dt, 11)
    t_eval_ndim = t_eval / sim.model.timescale.evaluate()

    # No external variables - Temperature solved as lumped model in pybamm
    # External variables could (and should) be used if battery thermal problem
    # Includes conduction with any other circuits or neighboring batteries
    # inp_and_ext.update(external_variables)
    inp_and_ext = inputs

    # Code to create mapped integrator
    integrator = solver.create_integrator(
        sim.built_model, inputs=inp_and_ext, t_eval=t_eval_ndim
    )
    integrator = integrator.map(Nspm, "thread", nproc)

    # Variables function for parallel evaluation
    casadi_objs = sim.built_model.export_casadi_objects(variable_names=variable_names)
    variables = casadi_objs["variables"]
    t, x, z, p = (
        casadi_objs["t"],
        casadi_objs["x"],
        casadi_objs["z"],
        casadi_objs["inputs"],
    )
    variables_stacked = casadi.vertcat(*variables.values())
    variables_fn = casadi.Function("variables", [t, x, z, p], [variables_stacked])
    variables_fn = variables_fn.map(Nspm, "thread", nproc)
    return integrator, variables_fn, t_eval


def solve(
    netlist=None,
    parameter_values=None,
    experiment=None,
    I_init=1.0,
    htc=None,
    initial_soc=0.5,
    nproc=12,
    output_variables=None,
):
    """
    Solves a pack simulation
    Parameters
    ----------
    netlist : pandas.DataFrame
        A netlist of circuit elements with format. desc, node1, node2, value.
        Produced by liionpack.read_netlist or liionpack.setup_circuit
    parameter_values : pybamm.ParameterValues class
        A dictionary of all the model parameters
    experiment : pybamm.Experiment class
        The experiment to be simulated. experiment.period is used to
        determine the length of each timestep.
    I_init : float, optional
        Initial guess for single battery current [A]. The default is 1.0.
    htc : float array, optional
        Heat transfer coefficient array of length Nspm. The default is None.
    initial_soc : float
        The initial state of charge for every battery. The default is 0.5
    nproc : int, optional
        Number of processes to start in parallel for mapping. The default is 12.
    output_variables : list, optional
        Variables to evaluate during solve. Must be a valid key in the
        model.variables
    Raises
    ------
    Exception
        DESCRIPTION.
    Returns
    -------
    output : ndarray shape [# variable, # steps, # batteries]
        simulation output array
    """

    if netlist is None or parameter_values is None or experiment is None:
        raise Exception("Please supply a netlist, paramater_values, and experiment")

    # Get netlist indices for resistors, voltage sources, current sources
    Ri_map = netlist["desc"].str.find("Ri") > -1
    V_map = netlist["desc"].str.find("V") > -1
    I_map = netlist["desc"].str.find("I") > -1
    Terminal_Node = np.array(netlist[I_map].node1)
    Nspm = np.sum(V_map)

    # Generate the protocol from the supplied experiment
    protocol = lp.generate_protocol_from_experiment(experiment)
    dt = experiment.period
    Nsteps = len(protocol)

    # Solve the circuit to initialise the electrochemical models
    V_node, I_batt = lp.solve_circuit(netlist)

    # Create battery simulation and update initial state of charge
    sim = lp.create_simulation(parameter_values, make_inputs=True)
    lp.update_init_conc(sim, SoC=initial_soc)

    # The simulation output variables calculated at each step for each battery
    # Must be a 0D variable i.e. battery wide volume average - or X-averaged for 1D model
    variable_names = [
        "Terminal voltage [V]",
        "Measured battery open circuit voltage [V]",
    ]
    if output_variables is not None:
        for out in output_variables:
            if out not in variable_names:
                variable_names.append(out)
        # variable_names = variable_names + output_variables
    Nvar = len(variable_names)

    # Storage variables for simulation data
    shm_i_app = np.zeros([Nsteps, Nspm], dtype=float)
    shm_Ri = np.zeros([Nsteps, Nspm], dtype=float)
    output = np.zeros([Nvar, Nsteps, Nspm], dtype=float)

    # Initialize currents in battery models
    shm_i_app[0, :] = I_batt * -1

    # Set up integrator
    integrator, variables_fn, t_eval = _create_casadi_objects(
        I_init, htc[0], sim, dt, Nspm, nproc, variable_names
    )

    # Step forward in time
    time = 0
    end_time = dt * Nsteps
    step_solutions = [None] * Nspm
    V_terminal = []
    record_times = []

    v_cut_lower = parameter_values["Lower voltage cut-off [V]"]
    v_cut_higher = parameter_values["Upper voltage cut-off [V]"]

    sim_start_time = ticker.time()

    for step in tqdm(range(Nsteps), desc='Solving Pack'):
        # Step the individual battery models
        step_solutions, var_eval = _mapped_step(
            sim.built_model,
            step_solutions,
            lp.build_inputs_dict(shm_i_app[step, :], htc),
            integrator,
            variables_fn,
            t_eval,
        )
        output[:, step, :] = var_eval

        time += dt

        # Calculate internal resistance and update netlist
        temp_v = output[0, step, :]
        temp_ocv = output[1, step, :]
        # temp_Ri = output[2, step, :]
        # This could be used instead of Equivalent ECM resistance which has
        # been changing definition
        temp_Ri = (temp_ocv - temp_v) / shm_i_app[step, :]
        # Make Ri more stable
        current_cutoff = np.abs(shm_i_app[step, :]) < 1e-6
        temp_Ri[current_cutoff] = 1e-12
        # temp_Ri = 1e-12
        shm_Ri[step, :] = temp_Ri

        netlist.loc[V_map, ("value")] = temp_ocv
        netlist.loc[Ri_map, ("value")] = temp_Ri
        netlist.loc[I_map, ("value")] = protocol[step]
        if time <= end_time:
            record_times.append(time)
            V_node, I_batt = lp.solve_circuit(netlist)
            V_terminal.append(V_node[Terminal_Node][0])
        if time < end_time:
            shm_i_app[step + 1, :] = I_batt[:] * -1
        # Stop if voltage limits are reached
        if np.any(temp_v < v_cut_lower):
            print("Low voltage limit reached")
            break
        if np.any(temp_v > v_cut_higher):
            print("High voltage limit reached")
            break

      

    # Collect outputs
    all_output = {}
    all_output["Time [s]"] = np.asarray(record_times)
    all_output["Pack current [A]"] = np.asarray(protocol[: step + 1])
    all_output["Pack terminal voltage [V]"] = np.asarray(V_terminal)
    all_output["Cell current [A]"] = shm_i_app[: step + 1, :]
    for j in range(Nvar):
        all_output[variable_names[j]] = output[j, : step + 1, :]

    toc = ticker.time()
    lp.logger.notice(
        "Solve circuit time " + str(np.around(toc - sim_start_time, 3)) + "s"
    )
    return all_output

Allow user to pass in simulation instead of parameter_values

This would probably help with issues #57 and #48 but might cause some confusion. The reason simulation was created behind the scenes was to setup the model with the thermal option and also to make the current and heat transfer coefficients inputs. Actually the thermal model doesn't need to be set by default and the inputs should be handled more generically allowing the user to specify anything as an input. The only think that needs to be an input is the current and even then that may be a sub-case for a more generic ensemble solve which we might want to allow for. For now liionpack does packs so this is probably a feature we can add (without too much trouble I expect)

Run Black formatting when PRs are merged

Make a drive cycle example

Different Initial condition for the cells in pack

I tried to populate different initial electrolyte concentration across the cells in the pack. I have an array of initial electrolyte concentration (initial_eleConc) similar to total heat transfer coefficient (htc). The changes I made to utils.py, simulations.py, solver_utils.py are given below.

Changed line 220, and added line 229 in utils.py.

# -*- coding: utf-8 -*-
"""
Created on Thu Sep 23 10:33:13 2021

@author: Tom
"""
from scipy.interpolate import interp1d, interp2d
import numpy as np
import pandas as pd
import os
import liionpack
from skspatial.objects import Plane
from skspatial.objects import Points

ROOT_DIR = os.path.dirname(os.path.abspath(liionpack.__file__))
CIRCUIT_DIR = os.path.join(ROOT_DIR, "circuits")
DATA_DIR = os.path.join(ROOT_DIR, "data")
INIT_FUNCS = os.path.join(ROOT_DIR, "init_funcs")


def interp_current(df):
    r"""
    Returns an interpolation function for current w.r.t time

    Parameters
    ----------
    df : pandas.DataFrame or Dict
        Contains data for 'Time' and 'Cells Total Current' from which to
        construct an interpolant function

    Returns
    -------
    f : function
        interpolant function of total cell current with time.

    """
    t = df["Time"]
    I = df["Cells Total Current"]
    f = interp1d(t, I)
    return f


def _z_from_plane(X, Y, plane):
    r"""
    Given X and Y and a plane provide Z
    X - temperature
    Y - flow rate
    Z - heat transfer coefficient

    Parameters
    ----------
    X : float (array)
        x-coordinate.
    Y : float (array)
        z-coordinate.
    plane : skspatial.object.Plane
        plane returned from read_cfd_data.

    Returns
    -------
    z : float (array)
        z-coordinate.

    """
    a, b, c = plane.normal
    d = plane.point.dot(plane.normal)
    z = (d - a * X - b * Y) / c
    return z


def _fit_plane(xv, yv, dbatt):
    r"""
    Private method to fit plane to CFD data

    Parameters
    ----------
    xv : ndarray
        temperature meshgrid points.
    yv : ndarray
        flow_rate meshgrid points.
    dbatt : ndarray
        cfd data for heat transfer coefficient.

    Returns
    -------
    plane : skspatial.object.Plane
        htc varies linearly with temperature and flow rate so relationship
        describes a plane

    """
    nx, ny = xv.shape
    pts = []
    for i in range(nx):
        for j in range(ny):
            pts.append([xv[i, j], yv[i, j], dbatt[i, j]])

    points = Points(pts)
    plane = Plane.best_fit(points, tol=1e-6)
    return plane


def read_cfd_data(data_dir=None, filename="cfd_data.xlsx", fit="linear"):
    r"""
    A very bespoke function to read heat transfer coefficients from an excel
    file

    Parameters
    ----------
    data_dir : str, optional
        Path to data file. The default is None. If unspecified the module
        liionpack.DATA_DIR folder will be used
    filename : str, optional
        The default is 'cfd_data.xlsx'.
    fit : str
        options are 'linear' (default) and 'interpolated'.
    Returns
    -------
    funcs : list
        an interpolant is returned for each cell in the excel file.

    """
    if data_dir is None:
        data_dir = liionpack.DATA_DIR
    fpath = os.path.join(data_dir, filename)
    ncells = 32
    flow_bps = np.array(pd.read_excel(fpath, sheet_name="massflow_bps", header=None))
    temp_bps = np.array(pd.read_excel(fpath, sheet_name="temperature_bps", header=None))
    xv, yv = np.meshgrid(temp_bps, flow_bps)
    data = np.zeros([len(temp_bps), len(flow_bps), ncells])
    fits = []
    for i in range(ncells):
        data[:, :, i] = np.array(
            pd.read_excel(fpath, sheet_name="cell" + str(i + 1), header=None)
        )
        # funcs.append(interp2d(xv, yv, data[:, :, i], kind='linear'))
        if fit == "linear":
            fits.append(_fit_plane(xv, yv, data[:, :, i]))
        elif fit == "interpolated":
            fits.append(interp2d(xv, yv, data[:, :, i], kind="linear"))

    return data, xv, yv, fits


def get_linear_htc(planes, T, Q):
    r"""
    A very bespoke function that is called in the solve process to update the
    heat transfer coefficients for every battery - assuming linear relation
    between temperature, flow rate and heat transfer coefficient.

    Parameters
    ----------
    planes : list
        each element of the list is a plane equation describing linear relation
        between temperature, flow rate and heat transfer coefficient.
    T : float array
        The temperature of each battery.
    Q : float
        The flow rate for the system.

    Returns
    -------
    htc : float
        Heat transfer coefficient for each battery.

    """
    ncell = len(T)
    htc = np.zeros(ncell)
    for i in range(ncell):
        htc[i] = _z_from_plane(T[i], Q, planes[i])
    return htc


def get_interpolated_htc(funcs, T, Q):
    r"""
    A very bespoke function that is called in the solve process to update the
    heat transfer coefficients for every battery

    Parameters
    ----------
    funcs : list
        each element of the list is an interpolant function.
    T : float array
        The temperature of each battery.
    Q : float
        The flow rate for the system.

    Returns
    -------
    htc : float
        Heat transfer coefficient for each battery.

    """
    ncell = len(T)
    htc = np.zeros(ncell)
    for i in range(ncell):
        htc[i] = funcs[i](T[i], Q)
    return htc


def build_inputs_dict(I_batt, htc,initial_eleConc):
    r"""
    Function to convert inputs and external_variable arrays to list of dicts
    As expected by the casadi solver. These are then converted back for mapped
    solving but stored individually on each returned solution.
    Can probably remove this process later

    Parameters
    ----------
    I_batt : float array
        The input current for each battery.
    htc : float array
        the heat transfer coefficient for each battery.

    Returns
    -------
    inputs_dict : list
        each element of the list is an inputs dictionary corresponding to each
        battery.


    """
    inputs_dict = []
    for i in range(len(I_batt)):
        inputs_dict.append(
            {
                # 'Volume-averaged cell temperature': T_batt[i],
                "Current": I_batt[i],
                "Total heat transfer coefficient [W.m-2.K-1]": htc[i],
                "Initial concentration in electrolyte [mol.m-3]": initial_eleConc[i],
            }
        )
    return inputs_dict

Added line 73 in simulations.py.

# -*- coding: utf-8 -*-
"""
Created on Wed Sep 22 15:37:51 2021

@author: tom
"""

import pybamm


def _current_function(t):
    r"""
    Internal function to make current an input parameter

    Parameters
    ----------
    t : float
        Dummy time parameter.

    Returns
    -------
    TYPE
        DESCRIPTION.

    """
    return pybamm.InputParameter("Current")


def create_simulation(parameter_values=None, experiment=None, make_inputs=False):
    r"""
    Create a PyBaMM simulation set up for interation with liionpack

    Parameters
    ----------
    parameter_values : pybamm.ParameterValues class
        DESCRIPTION. The default is None.
    experiment : pybamm.Experiment class
        DESCRIPTION. The default is None.
    make_inputs : bool, optional
        Changes "Current function [A]" and "Total heat transfer coefficient
        [W.m-2.K-1]" to be inputs that are controlled by liionpack.
        The default is False.

    Returns
    -------
    sim : pybamm.Simulation
        A simulation that can be solved individually or passed into the
        liionpack solve method

    """
    # Create the pybamm model
    model = pybamm.lithium_ion.SPMe(
        options={
            "thermal": "lumped",
        }
    )
    # geometry = model.default_geometry
    if parameter_values is None:
        # load parameter values and process model and geometry
        chemistry = pybamm.parameter_sets.Chen2020
        parameter_values = pybamm.ParameterValues(chemistry=chemistry)
    # Change the current function to be an input as this is set by the external circuit
    if make_inputs:
        parameter_values.update(
            {
                "Current function [A]": _current_function,
            }
        )
        parameter_values.update(
            {
                "Current": "[input]",
                "Total heat transfer coefficient [W.m-2.K-1]": "[input]",
                "Initial concentration in electrolyte [mol.m-3]":"[input]",
            },
            check_already_exists=False,
        )

    solver = pybamm.CasadiSolver(mode="safe")
    sim = pybamm.Simulation(
        model=model,
        experiment=experiment,
        parameter_values=parameter_values,
        solver=solver,
    )
    return sim


if __name__ == "__main__":
    sim = create_simulation()
    sim.solve([0, 1800])
    sim.plot()

Changed lines 86, 122,163,256,266 in solver_utils.py.

# -*- coding: utf-8 -*-
"""
Created on Thu Sep 23 10:44:31 2021

@author: Tom
"""
import casadi
import pybamm
import numpy as np
import time as ticker
import liionpack as lp


def _mapped_step(model, solutions, inputs_dict, integrator, variables, t_eval):
    r"""
    Internal function to process the model for one timestep in a mapped way.
    Mapped versions of the integrator and variables functions should already
    have been made.

    Parameters
    ----------
    model : pybamm.Model
        The built model
    solutions : list of pybamm.Solution objects for each battery
        Used to get the last state of the system and use as x0 and z0 for the
        casadi integrator
    inputs_dict : list of inputs_dict objects for each battery
        DESCRIPTION.
    integrator : mapped casadi.integrator
        Produced by _create_casadi_objects
    variables : mapped variables evaluator
        Produced by _create_casadi_objects
    t_eval : float array of times to evaluate
        Produced by _create_casadi_objects

    Returns
    -------
    sol : list
        solutions that have been stepped forward by one timestep
    var_eval : list
        evaluated variables for final state of system

    """
    len_rhs = model.concatenated_rhs.size
    N = len(solutions)
    if solutions[0] is None:
        # First pass
        x0 = casadi.horzcat(*[model.y0[:len_rhs] for i in range(N)])
        z0 = casadi.horzcat(*[model.y0[len_rhs:] for i in range(N)])
    else:
        x0 = casadi.horzcat(*[sol.y[:len_rhs, -1] for sol in solutions])
        z0 = casadi.horzcat(*[sol.y[len_rhs:, -1] for sol in solutions])
    # t_min = [0.0]*N
    t_min = 0.0
    inputs = []
    for temp in inputs_dict:
        inputs.append(casadi.vertcat(*[x for x in temp.values()] + [t_min]))
    ninputs = len(temp.values())
    inputs = casadi.horzcat(*inputs)
    # p = casadi.horzcat(*zip(inputs, external_variables, [t_min]*N))
    # inputs_with_tmin = casadi.vertcat(inputs, np.asarray(t_min))
    # Call the integrator once, with the grid
    timer = pybamm.Timer()
    tic = timer.time()
    casadi_sol = integrator(x0=x0, z0=z0, p=inputs)
    integration_time = timer.time()
    nt = len(t_eval)
    xf = casadi_sol["xf"]
    # zf = casadi_sol["zf"]
    sol = []
    xend = []
    for i in range(N):
        start = i * nt
        y_sol = xf[:, start:start + nt]
        xend.append(y_sol[:, -1])
        # Not sure how to index into zf - need an example
        sol.append(pybamm.Solution(t_eval, y_sol, model, inputs_dict[i]))
        sol[-1].integration_time = integration_time
    toc = timer.time()
    lp.logger.debug(f"Mapped step completed in {toc - tic}")
    xend = casadi.horzcat(*xend)
    var_eval = variables(0, xend, 0, inputs[0:ninputs, :])
    return sol, var_eval


def _create_casadi_objects(I_init, htc,initial_eleConc, sim, dt, Nspm, nproc, variable_names):
    r"""
    Internal function to produce the casadi objects in their mapped form for
    parallel evaluation

    Parameters
    ----------
    I_init : float
        initial guess for current of a battery (not used for simulation).
    htc : float
        initial guess for htc of a battery (not used for simulation).
    sim : pybamm.Simulation
        A PyBaMM simulation object that contains the model, parameter_values,
        solver, solution etc.
    dt : float
        The time interval for a single timestep. Fixed throughout the simulation
    Nspm : int
        Number of individual batteries in the pack.
    nproc : int
        Number of parallel processes to map to.
    variable_names : list
        Variables to evaluate during solve. Must be a valid key in the
        model.variables

    Returns
    -------
    integrator : mapped casadi.integrator
        Solves an initial value problem (IVP) coupled to a terminal value
        problem with differential equation given as an implicit ODE coupled
        to an algebraic equation and a set of quadratures
    variables_fn : mapped variables evaluator
        evaluates the simulation and output variables. see casadi function
    t_eval : float array of times to evaluate
        times to evaluate in a single step, starting at zero for each step

    """
    inputs = {"Current": I_init, "Total heat transfer coefficient [W.m-2.K-1]": htc,"Initial concentration in electrolyte [mol.m-3]": initial_eleConc}
    solver = sim.solver
    # solve model for 1 second to initialise the circuit
    t_eval = np.linspace(0, 1, 2)
    # Initial solution - this builds the model behind the scenes
    sim.solve(t_eval, inputs=inputs)
    # step model
    # Code to create mapped integrator
    t_eval = np.linspace(0, dt, 11)
    t_eval_ndim = t_eval / sim.model.timescale.evaluate()
    inp_and_ext = inputs
    # No external variables - Temperature solved as lumped model in pybamm
    # External variables could (and should) be used if battery thermal problem
    # Includes conduction with any other circuits or neighboring batteries
    # inp_and_ext.update(external_variables)

    integrator = solver.create_integrator(
        sim.built_model, inputs=inp_and_ext, t_eval=t_eval_ndim
    )
    integrator = integrator.map(Nspm, "thread", nproc)

    # Variables function for parallel evaluation
    casadi_objs = sim.built_model.export_casadi_objects(variable_names=variable_names)
    variables = casadi_objs["variables"]
    t, x, z, p = (
        casadi_objs["t"],
        casadi_objs["x"],
        casadi_objs["z"],
        casadi_objs["inputs"],
    )
    variables_stacked = casadi.vertcat(*variables.values())
    variables_fn = casadi.Function("variables", [t, x, z, p], [variables_stacked])
    variables_fn = variables_fn.map(Nspm, "thread", nproc)
    return integrator, variables_fn, t_eval


def solve(
    netlist=None,
    parameter_values=None,
    experiment=None,
    I_init=1.0,
    htc=None,initial_eleConc=None,
    initial_soc=0.5,
    nproc=12,
    output_variables=None,
):
    r"""
    Solves a pack simulation

    Parameters
    ----------
    netlist : pandas.DataFrame
        A netlist of circuit elements with format. desc, node1, node2, value.
        Produced by liionpack.read_netlist or liionpack.setup_circuit
    parameter_values : pybamm.ParameterValues class
        A dictionary of all the model parameters
    experiment : pybamm.Experiment class
        The experiment to be simulated. experiment.period is used to
        determine the length of each timestep.
    I_init : float, optional
        Initial guess for single battery current [A]. The default is 1.0.
    htc : float array, optional
        Heat transfer coefficient array of length Nspm. The default is None.
    initial_soc : float
        The initial state of charge for every battery. The default is 0.5
    nproc : int, optional
        Number of processes to start in parallel for mapping. The default is 12.
    output_variables : list, optional
        Variables to evaluate during solve. Must be a valid key in the
        model.variables

    Raises
    ------
    Exception
        DESCRIPTION.

    Returns
    -------
    output : ndarray shape [# variable, # steps, # batteries]
        simulation output array

    """

    if netlist is None or parameter_values is None or experiment is None:
        raise Exception("Please supply a netlist, paramater_values, and experiment")

    # Get netlist indices for resistors, voltage sources, current sources
    Ri_map = netlist["desc"].str.find("Ri") > -1
    V_map = netlist["desc"].str.find("V") > -1
    I_map = netlist["desc"].str.find("I") > -1

    Nspm = np.sum(V_map)

    protocol = lp.generate_protocol_from_experiment(experiment)
    dt = experiment.period
    Nsteps = len(protocol)

    # Solve the circuit to initialise the electrochemical models
    V_node, I_batt = lp.solve_circuit(netlist)

    sim = lp.create_simulation(parameter_values, make_inputs=True)
    lp.update_init_conc(sim, SoC=initial_soc)

    v_cut_lower = parameter_values["Lower voltage cut-off [V]"]
    v_cut_higher = parameter_values["Upper voltage cut-off [V]"]

    # The simulation output variables calculated at each step for each battery
    # Must be a 0D variable i.e. battery wide volume average - or X-averaged for 1D model
    variable_names = [
        "Terminal voltage [V]",
        "Measured battery open circuit voltage [V]",
        "Local ECM resistance [Ohm]",
    ]
    if output_variables is not None:
        for out in output_variables:
            if out not in variable_names:
                variable_names.append(out)
        # variable_names = variable_names + output_variables
    Nvar = len(variable_names)
    # Storage variables for simulation data
    shm_i_app = np.zeros([Nsteps, Nspm], dtype=float)
    shm_Ri = np.zeros([Nsteps, Nspm], dtype=float)
    output = np.zeros([Nvar, Nsteps, Nspm], dtype=float)

    # Initialize currents in battery models
    shm_i_app[0, :] = I_batt * -1

    time = 0
    # step = 0
    end_time = dt * Nsteps
    step_solutions = [None] * Nspm
    V_terminal = []
    record_times = []

    integrator, variables_fn, t_eval = _create_casadi_objects(
        I_init, htc[0],initial_eleConc[0], sim, dt, Nspm, nproc, variable_names
    )

    sim_start_time = ticker.time()

    for step in range(Nsteps):
        step_solutions, var_eval = _mapped_step(
            sim.built_model,
            step_solutions,
            lp.build_inputs_dict(shm_i_app[step, :], htc,initial_eleConc),
            integrator,
            variables_fn,
            t_eval,
        )
        output[:, step, :] = var_eval

        time += dt
        # Calculate internal resistance and update netlist
        temp_v = output[0, step, :]
        temp_ocv = output[1, step, :]
        temp_Ri = np.abs(output[2, step, :])
        shm_Ri[step, :] = temp_Ri

        netlist.loc[V_map, ("value")] = temp_ocv
        netlist.loc[Ri_map, ("value")] = temp_Ri
        netlist.loc[I_map, ("value")] = protocol[step]

        # print('Stepping time', np.around(ticker.time()-tic, 2), 's')
        if np.any(temp_v < v_cut_lower):
            print("Low V limit reached")
            break
        if np.any(temp_v > v_cut_higher):
            print("High V limit reached")
            break
        # step += 1
        if time <= end_time:
            record_times.append(time)
            V_node, I_batt = lp.solve_circuit(netlist)
            V_terminal.append(V_node.max())
        if time < end_time:
            shm_i_app[step + 1, :] = I_batt[:] * -1
    all_output = {}
    all_output["Time [s]"] = np.asarray(record_times)
    all_output["Pack current [A]"] = np.asarray(protocol[: step + 1])
    all_output["Pack terminal voltage [V]"] = np.asarray(V_terminal)
    all_output["Cell current [A]"] = shm_i_app[: step + 1, :]
    for j in range(Nvar):
        all_output[variable_names[j]] = output[j, : step + 1, :]

    toc = ticker.time()
    pybamm.logger.notice(
        "Solve circuit time " + str(np.around(toc - sim_start_time, 3)) + "s"
    )
    return all_output

The example code I ran is:

import liionpack as lp
import numpy as np
import pybamm

# Generate the netlist
netlist = lp.setup_circuit(Np=16, Ns=2, Rb=1e-4, Rc=1e-2, Ri=5e-2, V=3.2, I=80.0)

output_variables = [  
    'X-averaged total heating [W.m-3]',
    'Volume-averaged cell temperature [K]',
    'X-averaged negative particle surface concentration [mol.m-3]',
    'X-averaged positive particle surface concentration [mol.m-3]',
    'X-averaged electrolyte concentration [mol.m-3]',
    ]

# Heat transfer coefficients
htc = np.ones(32) * 10
#initial_eleConc=np.ones(32) * 1000.0
initial_eleConc=np.array([1000., 1000., 2000., 2000., 1000., 1000., 1000., 1000., 1000.,
       1000., 1000., 1000., 1000., 1000., 1000., 1000., 1000., 1000.,
       1000., 1000., 1000., 1000., 1000., 3000., 3000., 3000., 1000.,
       1000., 1000., 1000., 1000., 1000.])

# Cycling experiment, using PyBaMM
experiment = pybamm.Experiment(
    ["Charge at 50 A for 30 minutes", "Rest for 15 minutes", "Discharge at 50 A for 30 minutes", "Rest for 30 minutes"],
    period="10 seconds",
)

# PyBaMM parameters
chemistry = pybamm.parameter_sets.Chen2020
parameter_values = pybamm.ParameterValues(chemistry=chemistry)

# Solve pack
output = lp.solve(netlist=netlist,
                  parameter_values=parameter_values,
                  experiment=experiment,
                  output_variables=output_variables,
                  htc=htc,initial_eleConc=initial_eleConc)
lp.plot_output(output)

However, in the simulation the first value of initial electrolyte concentration, i.e., initial_eleConc[0] is populating across all the cells. Here is the output showing that

:
I think the issue is in line 256 of solver_utils.py, i.e.

integrator, variables_fn, t_eval = _create_casadi_objects(
        I_init, htc[0],initial_eleConc[0], sim, dt, Nspm, nproc, variable_names
    )

where it is only taking the first element of the array initial_eleConc[0]. I gave initial_eleConc array similar to htc(heat transfer coefficient) array. As htc array is going into lumped model, so every cell can have their own individual heat transfer coefficient. In contrast, the parameters in the base model seems not populating across the cells.

Set up benchmarks

Modular packs

Currently a pack is simply a single set of batteries in a combination of series and parallel with two busbars both with the terminal on the left hand side. This could represent a module in a larger pack and so needs to behave and interact as part of the larger system in a self contained and modular way.

Better plots

Can we implement something better for visualising large datasets

Finish Docs

move example notebooks to docs/examples and add to docs

I changed my mind about this - the examples look nice in the docs online and I don't think putting them inside the docs subdirectory hides them too badly

Include circuits, data, init_funcs in setup.py

Describe the bug

Installing liionpack using pip install -q git+https://github.com/pybamm-team/liionpack.git@main does not include circuits, data, init_funcs folders in the package

To Reproduce

Steps to reproduce the behaviour:
Run 02 Simulating a bespoke pack using an LTSpice netlist and custom heat transfer.ipynb from examples/notebooks

TO DO

Include these folders in setup.py

For Reference:
https://github.com/pybamm-team/PyBaMM/blob/5fbfec8a46945e99106cc166c283201513066425/setup.py#L136-L160

Set up google colab

Documenting multiple return values

Right now mkdocstrings doesn't support multiple return values or adding a name for the return values (hence the weird return values part in functions with > 1 return value).

This can be worked upon once the feature is added in.

Refer - mkdocstrings/mkdocstrings#301

ray multiprocessing test on github actions

See PR #87 for failing test

Getting this error running ray on github actions server. Could be something to do with how we're pip installing? @priyanshuone6 any ideas?
RuntimeError: The actor with name ray_actor failed to import on the worker. This may be because needed library dependencies are not installed in the worker environment

Small resting currents

The (calculated) internal resistance can increase massively around discontinuities and the starts of periods of rest. Investigate whether this is time-step dependent or whether Ri calculation might need to include a bound....

Test output of simulation by comparing to pybamm (integration test)

In some limit of an ideal pack (low resistances, high heat-transfer-coefficient?) the pack simulation should give an identical solution to a pybamm simulation of a single cell whose outputs have been scaled for the pack. This could be a good sanity check that the pack simulation is behaving as expected.

DFN model Run Time Error

I ran the example on liionpack by changing the Pybamm model from SPMe to DFN in 'simulations.py' file.
In 'simulations.py', I changed create_simulation function.
I replaced line #50: model = pybamm.lithium_ion.SPMe(options = {"thermal": "lumped",}) to
model = pybamm.lithium_ion.DFN(options = {"thermal": "lumped",})

I am getting a run time error as follows:

Newton/Linesearch algorithm failed to converge.
CasADi - 2021-10-19 09:06:28 WARNING("Exception raised: Error in Function::operator() for 'map3_F' [Map] at .../casadi/core/function.cpp:1368:
Error in Function::operator() for 'F' [IdasInterface] at .../casadi/core/function.cpp:1368:
.../casadi/interfaces/sundials/idas_interface.cpp:591: IDACalcIC returned "IDA_CONV_FAIL". Consult IDAS documentation.") [.../casadi/core/map.cpp:449]
Newton/Linesearch algorithm failed to converge.
CasADi - 2021-10-19 09:06:28 WARNING("Exception raised: Error in Function::operator() for 'map3_F' [Map] at .../casadi/core/function.cpp:1368:
Error in Function::operator() for 'F' [IdasInterface] at .../casadi/core/function.cpp:1368:
.../casadi/interfaces/sundials/idas_interface.cpp:591: IDACalcIC returned "IDA_CONV_FAIL". Consult IDAS documentation.") […/casadi/core/map.cpp:449]
RuntimeError: .../casadi/core/function_internal.hpp:1229: Evaluation failed

Solver warning for large pack configurations

I tested the example shown below with the following pack configurations:

Configuration	Elapsed time (min:sec)	Warnings
16p2s (32 cells)	0:08	none
32p10s (320 cells)	1:21	none
32p20s (640 cells)	3:56	none
16p40s (640 cells)	3:49	LinAlgWarning: Ill-conditioned matrix, X = solve(A, z).flatten()
16p96s (1536 cells)	18:43	LinAlgWarning: Ill-conditioned matrix, X = solve(A, z).flatten()

import liionpack as lp
import matplotlib.pyplot as plt
import numpy as np
import pybamm

# Parameters
Np = 16
Ns = 40

# Generate the netlist
netlist = lp.setup_circuit(Np=Np, Ns=Ns, Rb=1e-4, Rc=1e-2, Ri=5e-2, V=3.2, I=80.0)

output_variables = [
    'X-averaged total heating [W.m-3]',
    'Volume-averaged cell temperature [K]',
    'X-averaged negative particle surface concentration [mol.m-3]',
    'X-averaged positive particle surface concentration [mol.m-3]'
]

# Heat transfer coefficients
nbatt = Np * Ns
htc = np.ones(nbatt) * 10
print('nbatt =', nbatt)

# Cycling protocol
protocol = lp.generate_protocol()

# PyBaMM parameters
chemistry = pybamm.parameter_sets.Chen2020
parameter_values = pybamm.ParameterValues(chemistry=chemistry)

# Solve pack
output = lp.solve(netlist=netlist,
                  parameter_values=parameter_values,
                  protocol=protocol,
                  output_variables=output_variables,
                  htc=htc)

plt.show()

The ill conditioned warnings are related to the scipy.linalg.solve function in netlist_utils.py which solves the electrical circuit for node voltage and battery current. The A matrix used by the solver for X = solve(A, z).flatten() appears to be a banded matrix. Therefore, it may be appropriate to use the scipy.linalg.solve_banded function to solve the system of equations representing the circuit. This may also improve the solution time of the pack simulation. But I'm not sure if all pack configurations would generate a banded A matrix for the circuit equations. The code may need to check if the configuration is suitable for the scipy.linalg.solve or for the scipy.linalg.solve_banded function.

Separate plot from solve

Separate plotting features from solve so manipulation of the plots can happen without solving the model every time.

Prep

Read all this and set-up rest of todo
https://joss.readthedocs.io/en/latest/submitting.html

Turn off latex printing in Ipython console after drawing circuit

Adding external resistance option

For cell modelling it would be helpful to add a term that takes into account external series resistances such as tab weld resistances.
Would it be possible to add an optional resistance term to the simulation?

Thanks in advance!

Experiments instead of Protocols

Protocols is very quick and dirty just to get some simple cycling behaviour working. PyBaMM has much nicer Experiment class that hopefully we can leverage here. Although it may be simpler to cut and paste I would like to avoid that

module 'skspatial.objects' not found

Hi @TomTranter, you use a module named 'skspatial.objects' in the liionpack package - where can I get this from? Is it linked to scipy? I can't seem to find this module anywhere...
Many thanks in advance!

using parallel jax solver

@TomTranter : here is an example script that runs 20 instances of the spm model in parallel using jax. Let me know if this is what you need.

import time
import pybamm
import numpy as np
import jax
import matplotlib.pylab as plt

import os
# specify 20 logical devices for execution
ncpu = 20
os.environ['XLA_FLAGS'] = (
    '--xla_force_host_platform_device_count={}'.format(ncpu)
)

# print out the available devices

print('devices', jax.devices())

pybamm.set_logging_level("INFO")
model = pybamm.lithium_ion.SPM()
model.convert_to_format = "jax"
model.events = []

# create geometry
geometry = model.default_geometry

# load parameter values and process model and geometry
param = model.default_parameter_values
parameter = "Electrode height [m]"
value = param[parameter]
param.update({parameter: "[input]"})
param.process_model(model)
param.process_geometry(geometry)

# set mesh
mesh = pybamm.Mesh(geometry, model.default_submesh_types, model.default_var_pts)

# discretise model
disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
disc.process_model(model)

# solve model for 1 hour
t_eval = np.linspace(0, 3600, 100)
solver = pybamm.JaxSolver()

# the model is setup and compiled (expensive) during this call
solution = solver.solve(
    model, t_eval,
    inputs={parameter: value},
)

# create a new jax function that can take an array of inputs
mapped_jax_function = jax.pmap(
    solver.get_solve(model, t_eval),
)

# now our input is an array of "parameter",
# size needs to be the same as the number
# of devices
inputs_array = {
    parameter: jax.numpy.linspace(value / 10, value * 10, ncpu)
}

# the multiple inputs are executed in parallel here,
# the result is a 3d array of shape (Ni, Ns, Nt), where
# Ni is the size of the parameter array, Ns is the size of the
# full state vector, and Nt is the number of timesteps in t_eval
print('running in parallel')
tic = time.perf_counter()
result_array = mapped_jax_function(inputs_array)

# access the result so its actually computed
print(result_array[0, 0, 0])
toc = time.perf_counter()
print('time elapsed: {} sec', toc - tic)

GitHub Action workflow to check conda environment

Can someone add a GitHub Action to build and check the conda environment? I don't have permission to update the workflow in this repository but I believe it would be something like this:

  build:
    needs: style
    runs-on: ubuntu-latest

    steps:
    - uses: actions/checkout@v2
    - name: Set up Python 3.9
      uses: actions/setup-python@v2
      with:
        python-version: 3.9
    - name: Install dependencies
      run: |
        $CONDA/bin/conda env update --file environment.yml --name base

Add more debug info to solve

Liionpack uses the opposite sign convention from pybamm

In pybamm, discharge current is positive and charge current is negative.
In liionpack, it is the other way round. This might cause confusion. Can I just fix as part of #9?