Obtain solution results without using ProcessedVariable

[wigging]

The Solution object has a y property that represents the values of the solution. But how do I use this property to obtain values such as voltage, current, time, etc.? For example, the voltage values can be obtained from the solution's ProcessedVariable using solution['Terminal voltage [V]'].entries. How would I get the voltage results without using the processed variable object?

The example below calculates the solution at each time step then stores the time, current, and voltage. Using cProfile, the processed variables account for 25% of the time it takes to execute the script. It would be less expensive if there is a way to get the results out of the solution without using the ProcessedVariable objects at each time step.

import numpy as np
import pybamm
import time

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

    # Parameters
    discharge = np.linspace(1, 4, num=10)
    dt = 1
    tf = 4000

    # Load model and create geometry
    model = pybamm.lithium_ion.SPMe()
    geometry = model.default_geometry

    # Load parameter values then process model and geometry
    param = model.default_parameter_values
    param['Current function [A]'] = '[input]'
    param.process_model(model)
    param.process_geometry(geometry)

    # Set the 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)

    # Solver for model
    solver = pybamm.CasadiSolver()

    # Solve at each time step
    times = []
    currents = []
    voltages = []

    for dis in discharge:
        t = 0
        sol_t = []
        sol_c = []
        sol_v = []

        t_eval = np.linspace(0, 1, 2)
        sol = solver.solve(model, t_eval, inputs={"Current function [A]": dis})
        sol_t.append(sol['Time [s]'].entries[0])
        sol_c.append(sol['Current [A]'].entries[0])
        sol_v.append(sol['Terminal voltage [V]'].entries[0])

        while t < tf:
            sol = solver.step(sol, model, dt=dt, npts=2, inputs={"Current function [A]": dis}, save=False)
            sol_t.append(sol['Time [s]'].entries[0])
            sol_c.append(sol['Current [A]'].entries[0])
            sol_v.append(sol['Terminal voltage [V]'].entries[0])
            t += dt

        times.append(sol_t)
        currents.append(sol_c)
        voltages.append(sol_v)

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

if __name__ == '__main__':
    main()

[valentinsulzer]

At best you could look inside the ProcessedVariable class and do that manually. But what's the use case?

[wigging]

When running lots of simulations in parallel and getting results at each time step, I noticed there is some overhead associated with the ProcessedVariable. I was wondering if there is a faster way to get the results instead of calling solution['Terminal voltage [V]'].entries at each time step. I thought the y property from the solution object might provide a quicker way to get the results.

[valentinsulzer]

There are (at least) two overheads:

1. Converting the variable to casadi before post-processing

   PyBaMM/pybamm/solvers/solution.py

   Lines 456 to 481 in dd77bec

   	for model, ys, inputs, var_pybamm in zip(
   	self.all_models, self.all_ys, self.all_inputs, vars_pybamm
   	):
   	if key in model._variables_casadi:
   	var_casadi = model._variables_casadi[key]
   	else:
   	t_MX = casadi.MX.sym("t")
   	y_MX = casadi.MX.sym("y", ys.shape[0])
   	symbolic_inputs_dict = {
   	key: casadi.MX.sym("input", value.shape[0])
   	for key, value in inputs.items()
   	}
   	symbolic_inputs = casadi.vertcat(
   	*[p for p in symbolic_inputs_dict.values()]
   	)
   	# Convert variable to casadi
   	# Make all inputs symbolic first for converting to casadi
   	var_sym = var_pybamm.to_casadi(
   	t_MX, y_MX, inputs=symbolic_inputs_dict
   	)
   	var_casadi = casadi.Function(
   	"variable", [t_MX, y_MX, symbolic_inputs], [var_sym]
   	)
   	model._variables_casadi[key] = var_casadi

   
   The way to avoid this overhead is to make sure the same model object is passed to solution each time (so that the model only gets set up once)
2. Computing the solution at every time step, which is currently done in a loop but could be parallelized

   PyBaMM/pybamm/solvers/processed_variable.py

   Lines 112 to 119 in dd77bec

   	for ts, ys, inputs, base_var_casadi in zip(
   	self.all_ts, self.all_ys, self.all_inputs_casadi, self.base_variables_casadi
   	):
   	for inner_idx, t in enumerate(ts):
   	t = ts[inner_idx]
   	y = ys[:, inner_idx]
   	entries[idx] = base_var_casadi(t, y, inputs).full()[0, 0]
   	idx += 1

   
   If you only need the solution at the first or last step, you could use solution.first_state or solution.first_state

If you want to write your own version, all the code for post-processing is in Solution and ProcessedVariable

[wigging]

I still don't understand how this translates to a voltage or current value. What are all the casadi variables? Are they objects generated and used by the Casadi solver?

[wigging]

I updated my question with some example code. Based on the example, how do I use the Casadi objects to get the time step results without relying on the processed variable objects?

[valentinsulzer]

This is all post-processing after the casadi solver has finished. It evaluates the complicated expressions for e.g. voltage, which are not explicit states of the model

We use casadi because it is faster for post-processing. All of the code you need is in Solution.update, ProcessedVariable.__init__, and ProcessedVariable.initialise_0D (each called once in that order).

[valentinsulzer]

y contains variables of the model like concentrations and potentials, then there is a post-processing step to get the current and voltage out. Time is solution.t