I’m relative new to Pybamm. I’m working on a parallel-in-time methods.
I’ve noticed some differences in the starting times in the solution of an experiment with more than 1 cycle. Here is a small example:
# -*- coding: utf-8 -*-
"""
Spyder Editor
This is a temporary script file.
"""
import pybamm
import numpy as np
nCycles = 10
model = pybamm.lithium_ion.DFN({"working electrode":"positive"})
var_pts = {
"x_n": 1, # points in the negative electrode
"x_s": 30, # points in the separator
"x_p": 30, # points in the positive electrode
"r_p": 100 # points in the radius of positive electrode
}
parameter_values=pybamm.ParameterValues("Xu2019")
experiment = pybamm.Experiment(
[
(
#"Rest for 5 minutes",
"Discharge at 1C until 3. V",
"Charge at 1C until 4.2 V",
)
]
* nCycles
)
solver = pybamm.IDAKLUSolver()
sim_seq = pybamm.Simulation(
model,
experiment=experiment,
parameter_values=parameter_values,
solver=solver,
var_pts=var_pts,)
sol_seq = sim_seq.solve()
cycle = 1
for sol_before, sol_after in zip(sol_seq.cycles[:-1], sol_seq.cycles[1:]):
print(f"Cycle = {cycle}\n")
print("Previous solution ends at t: ", sol_before.t[-1]," == last_state.t", sol_before.last_state.t[0] == sol_before.t[-1])
print("Subsequent solution starts at t:", sol_after.t[0]," == first_state.t", sol_after.first_state.t[0]==sol_after.t[0])
print(f"Differences in initializing: {np.abs(sol_before.t[-1] - sol_after.t[0])}")
y_before = np.asarray(sol_before.last_state.y).reshape(-1)
y_after = np.asarray(sol_after.first_state.y).reshape(-1)
err_y = np.max(np.abs(y_before - y_after))
print(f"State differeneces: {err_y}\n")
cycle+=1
tsEnd = [sol.last_state.t for sol in sol_seq.cycles]
tsStart = [sol.first_state.t for sol in sol_seq.cycles]
print(np.squeeze(tsEnd[:-1])-np.squeeze(tsStart[1:]))
Hi Geri, this is intentional, the relevent lines are solvers/base_solver.py:1446-1456. We find the next double point number ofter the end of the previous time step so that there are no repeated time points in the solution object. I would imagine this is done so that the interpolation between timepoints does not break, but am not sure without more digging. Why do you want the starting times to line up exactly?
Initially I wanted to test the results between sequential cycling (all cycles at once) and stepwise (one cycle after another). Then i noticed that the values of the state variables and starting times are not matching.
I just wanted to start somewhere to build a propagator from one state to another. That propagator has to match the sequential solver if they are of the same “setup”.
since you are interested in parallel in time methods: I have recently implemented batching in diffsol, an ODE solver I’m writing (GitHub - martinjrobins/diffsol: ODE solver library in Rust · GitHub). I think batching would be a nice way to implement POT methods, but have not had a chance to try it out yet.
Thank you!
Maybe you can confirm then: this correction is then done at each Experiment step ?
I guess the experiment*number_of_cycles in Pybamm is just “one cycle” and the correction is done at each step ?
yup, each experimental step is a separate solve and the solutions of each subsequent solve need to be concatenated together, so the correction is done between each step
Hello!
I have another question to the same implementation:
experiment_step = pybamm.Experiment(
[
(
#"Rest for 5 minutes",
"Discharge at 1C until 3. V",
"Charge at 1C until 4.2 V",
)
]
* 1
)
sim_step = pybamm.Simulation(
model,
experiment=experiment_step,
parameter_values=parameter_values,
solver=solver,
var_pts=var_pts,)
# solve stepwise
init_sol = None
solution_step = []
for cycle in range(nCycles):
sol = sim_step.solve(starting_solution=init_sol)
cycle_sol = sol.cycles[-1]
solution_step.append(cycle_sol)
init_sol = cycle_sol.last_state
print("Differences in states between step and sequential solutions for the first 2 cycles\n")
sol_step_0 = solution_step[0]
sol_step_1 = solution_step[1]
sol_seq_0 = sol_seq.cycles[0]
sol_seq_1 = sol_seq.cycles[1]
n_diff = sol_seq_0.all_models[0].concatenated_rhs.size
y_size = sol_seq_0.first_state.y.shape[0]
n_alg = sol_seq_0.all_models[0].concatenated_algebraic.size
print("Delta t_minus:", sol_step_0.last_state.t - sol_seq_0.last_state.t)
print("Delta t_plus:", sol_step_1.first_state.t - sol_seq_1.first_state.t)
print("Delta y_minus:", np.linalg.norm(sol_step_0.last_state.y[:,-1] - sol_seq_0.last_state.y[:,-1], ord = np.inf))
print("Delta y_plus:", np.linalg.norm(sol_step_1.first_state.y[:,-1] - sol_seq_1.first_state.y[:,-1], ord = np.inf))
print("Delta y[DIFFERENTIAL PART]_plus:", np.linalg.norm(sol_step_1.first_state.y[:n_diff,-1] - sol_seq_1.first_state.y[:n_diff,-1], ord = np.inf))
print("Delta y[ALGEBRAIC PART]_plus:", np.linalg.norm(sol_step_1.first_state.y[n_diff:,-1] - sol_seq_1.first_state.y[n_diff:,-1], ord = np.inf))
OUTPUT:
Differences in states between step and sequential solutions for the first 2 cycles
Delta t_minus: [0.]
Delta t_plus: [0.]
Delta y_minus: 0.0
Delta y_plus: 3.193136033363686e-11
Delta y[DIFFERENTIAL PART]_plus: 2.4807933485249123e-13
Delta y[ALGEBRAIC PART]_plus: 3.193136033363686e-11
I don’t see where I might have done some calculations wrong… Under Assumption that it is correct ( haha ) I expected the two computations deliver same results. Could someone please help me ?
In both cycles the time correction is the same, I expect the starting condition for the second cycle to be identical for both, but it is not…
(*)_minus and (*)_plus mean left and right limits of the (*)