Skip to content

Commit 16468ae

Browse files
rtimmsrtimms
committed
#853 fix event eval to nan
1 parent f1af7a0 commit 16468ae

File tree

3 files changed

+196
-15
lines changed

3 files changed

+196
-15
lines changed

‎pouch_cell.py

+163
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,163 @@
1+
import pybamm
2+
import numpy as np
3+
import matplotlib.pyplot as plt
4+
5+
pybamm.set_logging_level("INFO")
6+
7+
# load model
8+
options = {
9+
"current collector": "potential pair",
10+
"dimensionality": 2,
11+
"thermal": "x-lumped",
12+
# The below option replaces the PDEs in the particles with ODEs under the
13+
# assumption of fast diffusion (so that the concentration is uniform within
14+
# each particle). This will speed up the simulation and is an OK assumption
15+
# for a lot of cases. Uncomment it to switch it on/off.
16+
"particle": "fast diffusion",
17+
}
18+
model = pybamm.lithium_ion.SPM(options)
19+
20+
# parameters can be updated here
21+
param = model.default_parameter_values
22+
23+
# set mesh points
24+
var = pybamm.standard_spatial_vars
25+
var_pts = {
26+
var.x_n: 10,
27+
var.x_s: 10,
28+
var.x_p: 10,
29+
var.r_n: 10,
30+
var.r_p: 10,
31+
var.y: 10,
32+
var.z: 10,
33+
}
34+
35+
# solver
36+
# casadi fast mode is pretty quick, but doesn't support events out of the box
37+
solver = pybamm.CasadiSolver(atol=1e-6, rtol=1e-3, root_tol=1e-3, mode="fast")
38+
# KLU is sometimes better for bigger problems and supports events
39+
# solver = pybamm.IDAKLUSolver(atol=1e-6, rtol=1e-3, root_tol=1e-3, root_method="hybr")
40+
# KLU performs better if you convert the final model to a python function
41+
if isinstance(solver, pybamm.IDAKLUSolver):
42+
model.convert_to_format = "python"
43+
44+
# simulation object
45+
simulation = pybamm.Simulation(
46+
model, parameter_values=param, var_pts=var_pts, solver=solver
47+
)
48+
49+
# build simulation
50+
# by default pybamm performs some checks on the discretised model. this can
51+
# be a little slow for bigger problems, so you can turn if off. if you start to
52+
# get obscure errors when you change things, set this to True and the error should
53+
# get caught sooner and give a more informative message
54+
simulation.build(check_model=False)
55+
56+
# solve simulation
57+
t_eval = np.linspace(0, 3600, 100) # time in seconds
58+
simulation.solve(t_eval=t_eval)
59+
solution = simulation.solution
60+
61+
# plotting
62+
# TO DO: 2+1D automated plotting
63+
64+
# post-process variables
65+
phi_s_cn = solution["Negative current collector potential [V]"]
66+
phi_s_cp = solution["Positive current collector potential [V]"]
67+
I = solution["Current collector current density [A.m-2]"]
68+
T = solution["X-averaged cell temperature [K]"]
69+
70+
# get y and z points for plotting (these are non-dimensional)
71+
l_y = phi_s_cp.y_sol[-1]
72+
l_z = phi_s_cp.z_sol[-1]
73+
y_plot = np.linspace(0, l_y, 21)
74+
z_plot = np.linspace(0, l_z, 21)
75+
76+
# Can multiply by L_z to get dimensional y and z. Note that both y and z are
77+
# scaled with L_z
78+
L_z = param.evaluate(pybamm.standard_parameters_lithium_ion.L_z)
79+
y_plot_dim = np.linspace(0, l_y, 21) * L_z
80+
z_plot_dim = np.linspace(0, l_z, 21) * L_z
81+
82+
83+
# define plotting function
84+
def plot(time_in_seconds):
85+
fig, ax = plt.subplots(figsize=(15, 8))
86+
plt.tight_layout()
87+
plt.subplots_adjust(left=-0.1)
88+
89+
# get non-dim time
90+
tau = param.evaluate(pybamm.standard_parameters_lithium_ion.tau_discharge)
91+
t_non_dim = time_in_seconds / tau
92+
93+
# negative current collector potential
94+
plt.subplot(221)
95+
phi_s_cn_plot = plt.pcolormesh(
96+
y_plot_dim,
97+
z_plot_dim,
98+
np.transpose(
99+
phi_s_cn(y=y_plot, z=z_plot, t=t_non_dim)
100+
), # accepts non-dim values
101+
shading="gouraud",
102+
)
103+
plt.axis([0, l_y * L_z, 0, l_z * L_z])
104+
plt.xlabel(r"$y$ [m]")
105+
plt.ylabel(r"$z$ [m]")
106+
plt.title(r"$\phi_{s,cn}$ [V]")
107+
plt.set_cmap("cividis")
108+
plt.colorbar(phi_s_cn_plot)
109+
110+
# positive current collector potential
111+
plt.subplot(222)
112+
phi_s_cp_plot = plt.pcolormesh(
113+
y_plot_dim,
114+
z_plot_dim,
115+
np.transpose(
116+
phi_s_cp(y=y_plot, z=z_plot, t=t_non_dim)
117+
), # accepts non-dim values
118+
shading="gouraud",
119+
)
120+
plt.axis([0, l_y * L_z, 0, l_z * L_z])
121+
plt.xlabel(r"$y$ [m]")
122+
plt.ylabel(r"$z$ [m]")
123+
plt.title(r"$\phi_{s,cp}$ [V]")
124+
plt.set_cmap("viridis")
125+
plt.colorbar(phi_s_cp_plot)
126+
127+
# through-cell current
128+
plt.subplot(223)
129+
I_plot = plt.pcolormesh(
130+
y_plot_dim,
131+
z_plot_dim,
132+
np.transpose(I(y=y_plot, z=z_plot, t=t_non_dim)), # accepts non-dim values
133+
shading="gouraud",
134+
)
135+
plt.axis([0, l_y * L_z, 0, l_z * L_z])
136+
plt.xlabel(r"$y$ [m]")
137+
plt.ylabel(r"$z$ [m]")
138+
plt.title(r"$I$ [A.m-2]")
139+
plt.set_cmap("plasma")
140+
plt.colorbar(I_plot)
141+
# temperature
142+
plt.subplot(224)
143+
T_plot = plt.pcolormesh(
144+
y_plot_dim,
145+
z_plot_dim,
146+
np.transpose(T(y=y_plot, z=z_plot, t=t_non_dim)), # accepts non-dim values
147+
shading="gouraud",
148+
)
149+
plt.axis([0, l_y * L_z, 0, l_z * L_z])
150+
plt.xlabel(r"$y$ [m]")
151+
plt.ylabel(r"$z$ [m]")
152+
plt.title(r"$T$ [K]")
153+
plt.set_cmap("inferno")
154+
plt.colorbar(T_plot)
155+
156+
plt.subplots_adjust(
157+
top=0.92, bottom=0.15, left=0.10, right=0.9, hspace=0.5, wspace=0.5
158+
)
159+
plt.show()
160+
161+
162+
# call plot
163+
plot(1000)

‎pybamm/models/submodels/current_collector/potential_pair.py

+1-2
Original file line numberDiff line numberDiff line change
@@ -69,13 +69,12 @@ def set_algebraic(self, variables):
6969
def set_initial_conditions(self, variables):
7070

7171
applied_current = self.param.current_with_time
72-
cc_area = self._get_effective_current_collector_area()
7372
phi_s_cn = variables["Negative current collector potential"]
7473
i_boundary_cc = variables["Current collector current density"]
7574

7675
self.initial_conditions = {
7776
phi_s_cn: pybamm.Scalar(0),
78-
i_boundary_cc: applied_current / cc_area,
77+
i_boundary_cc: applied_current,
7978
}
8079

8180

‎pybamm/solvers/casadi_solver.py

+32-13
Original file line numberDiff line numberDiff line change
@@ -145,13 +145,13 @@ def _integrate(self, model, t_eval, inputs=None):
145145
# Try to integrate in global steps of size dt_max (arbitrarily taken
146146
# to be half of the range of t_eval, up to a maximum of 1 hour)
147147
t_f = t_eval[-1]
148-
# dt_max must be at least as big as the the biggest step in t_eval
149-
# to avoid an empty integration window below
150-
dt_eval_max = np.max(np.diff(t_eval))
151-
dt_max = np.max([(t_f - t) / 2, dt_eval_max])
152-
# maximum dt_max corresponds to 1 hour
153148
t_one_hour = 3600 / model.timescale_eval
154-
dt_max = np.min([dt_max, t_one_hour])
149+
dt_max = np.min([(t_f - t) / 2, t_one_hour])
150+
# dt_max must be at least as big as the the biggest step in t_eval
151+
# (multiplied by some tolerance, here 0.01) to avoid an empty
152+
# integration window below
153+
dt_eval_max = np.max(np.diff(t_eval)) * 1.01
154+
dt_max = np.max([dt_max, dt_eval_max])
155155
while t < t_f:
156156
# Step
157157
solved = False
@@ -167,7 +167,7 @@ def _integrate(self, model, t_eval, inputs=None):
167167
# Try to solve with the current global step, if it fails then
168168
# halve the step size and try again.
169169
try:
170-
# When not in DEBUG mode (level=10), supress the output
170+
# When not in DEBUG mode (level=10), suppress the output
171171
# from the failed steps in the solver
172172
if (
173173
pybamm.logger.getEffectiveLevel() == 10
@@ -221,13 +221,32 @@ def _integrate(self, model, t_eval, inputs=None):
221221
# loop over events to compute the time at which they were triggered
222222
t_events = [None] * len(active_events)
223223
for i, event in enumerate(active_events):
224-
t_events[i] = brentq(
225-
lambda t: event(t, y_sol(t), inputs),
226-
current_step_sol.t[0],
227-
current_step_sol.t[-1],
228-
)
224+
225+
def event_fun(t):
226+
return event(t, y_sol(t), inputs)
227+
228+
if np.isnan(event_fun(current_step_sol.t[-1])[0]):
229+
# bracketed search fails if f(a) or f(b) is NaN, so we
230+
# need to find the times for which we can evaluate the event
231+
times = [
232+
t
233+
for t in current_step_sol.t
234+
if event_fun(t)[0] == event_fun(t)[0]
235+
]
236+
else:
237+
times = current_step_sol.t
238+
# skip if sign hasn't changed
239+
if np.sign(event_fun(times[0])) != np.sign(
240+
event_fun(times[-1])
241+
):
242+
t_events[i] = brentq(
243+
lambda t: event_fun(t), times[0], times[-1]
244+
)
245+
else:
246+
t_events[i] = np.nan
247+
229248
# t_event is the earliest event triggered
230-
t_event = np.min(t_events)
249+
t_event = np.nanmin(t_events)
231250
y_event = y_sol(t_event)
232251

233252
# return truncated solution

0 commit comments

Comments
 (0)