Skip to content

Commit

Permalink
allow for multi-direction experiment protocols
Browse files Browse the repository at this point in the history
  • Loading branch information
@mleot
mleot committed Sep 17, 2024
1 parent 6681c42 commit 70a647e
Show file tree
Hide file tree
Showing 3 changed files with 874 additions and 84 deletions.
675 changes: 675 additions & 0 deletions docs/examples/04 Running a multi step simulation.ipynb
View file

Large diffs are not rendered by default.

179 changes: 98 additions & 81 deletions jellybamm/__liionsolve__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,9 +59,9 @@ def run_simulation_lp(parameter_values, experiment, initial_soc, project):
# sorted_res_Ts = net["throat.spm_resistor_order"][res_Ts].argsort()
# print("Total Electrode Height", np.around(np.sum(electrode_heights), 2), "m")
# Take I_app from first command of the experiment
protos, terms, steps = lp.generate_protocol_from_experiment(experiment)
if len(protos) > 1:
raise ValueError("Experiment must be a single step currently")
protos, terms, types = lp.generate_protocol_from_experiment(experiment)
# if len(protos) > 1:
# raise ValueError("Experiment must be a single step currently")
I_app = protos[0][0]
I_typical = I_app / Nspm

Expand Down Expand Up @@ -130,75 +130,82 @@ def run_simulation_lp(parameter_values, experiment, initial_soc, project):
}
###########################################################################
# Initialisation
experiment_init = pybamm.Experiment(
[
f"Discharge at {I_app} A for 4 seconds or until 0V",
],
period="1 second",
)
proto_init, term_init, type_init = lp.generate_protocol_from_experiment(experiment_init)
# Solve the pack
manager = lp.CasadiManager()
manager.solve(
netlist=netlist,
sim_func=lp.thermal_external,
parameter_values=parameter_values,
experiment=experiment_init,
output_variables=output_variables,
inputs=inputs,
nproc=max_workers,
initial_soc=initial_soc,
setup_only=True,
)
Qvar = "Volume-averaged total heating [W.m-3]"
Qid = np.argwhere(np.asarray(manager.variable_names) == Qvar).flatten()[0]
lp.logger.notice("Starting initial step solve")
vlims_ok = True
tic = ticker.time()
netlist["power_loss"] = 0.0
# plt.figure()
skip_vcheck = True
manager.global_step = 0
with tqdm(total=len(proto_init[0]), desc="Initialising simulation") as pbar:
step = 0
while step < len(proto_init[0]):
###################################################################
updated_inputs = {"Input temperature [K]": spm_temperature}
vlims_ok = manager._step(step, proto_init[0], term_init[0], type_init[0], updated_inputs, skip_vcheck)
###################################################################
# Apply Heat Sources
Q_tot = manager.output[Qid, step, :]
Q = get_cc_power_loss(net, netlist)
# To do - Get cc heat from netlist
# Q_ohm_cc = net.interpolate_data("pore.cc_power_loss")[res_Ts]
# Q_ohm_cc /= net["throat.volume"][res_Ts]
# key = "Volume-averaged Ohmic heating CC [W.m-3]"
# vh[key][outer_step, :] = Q_ohm_cc[sorted_res_Ts]
Q[res_Ts] += Q_tot
jellybamm.apply_heat_source_lp(project, Q)
# Calculate Global Temperature
jellybamm.run_step_transient(
project, dim_time_step, T0, cp, rho, thermal_third
)
# Interpolate the node temperatures for the SPMs
spm_temperature = phase.interpolate_data("pore.temperature")[res_Ts]
# T_non_dim_spm = fT_non_dim(parameter_values, spm_temperature)
###################################################################
if vlims_ok:
step += 1
pbar.update(1)
temp_Ri = np.array(netlist.loc[manager.Ri_map].value)
# plt.scatter(np.arange(len(temp_Ri)), temp_Ri, label=str(step))
else:
break
# plt.legend()
manager.step = step
toc = ticker.time()
lp.logger.notice("Initial step solve finished")
lp.logger.notice("Total stepping time " + str(np.around(toc - tic, 3)) + "s")
lp.logger.notice(
"Time per step " + str(np.around((toc - tic) / manager.Nsteps, 3)) + "s"
)
def initialize_simulation(I_app, netlist, project, phase, spm_temperature):
if I_app < 0:
experiment_string = f"Discharge at {abs(I_app)} A for 4 seconds or until 0V"
elif I_app > 0:
experiment_string = f"Charge at {I_app} A for 4 seconds or until 4.5V"
else:
raise ValueError("I_app must be non-zero")
experiment_init = pybamm.Experiment([experiment_string], period="1 second")

proto_init, term_init, type_init = lp.generate_protocol_from_experiment(experiment_init)
# Solve the pack
tmp_manager = lp.CasadiManager()
tmp_manager.solve(
netlist=netlist,
sim_func=lp.thermal_external,
parameter_values=parameter_values,
experiment=experiment_init,
output_variables=output_variables,
inputs=inputs,
nproc=max_workers,
initial_soc=initial_soc,
setup_only=True,
)
Qvar = "Volume-averaged total heating [W.m-3]"
Qid = np.argwhere(np.asarray(tmp_manager.variable_names) == Qvar).flatten()[0]
lp.logger.notice("Starting initial step solve")
vlims_ok = True
tic = ticker.time()
netlist["power_loss"] = 0.0
# plt.figure()
skip_vcheck = True
tmp_manager.global_step = 0
with tqdm(total=len(proto_init[0]), desc="Initialising simulation") as pbar:
step = 0
while step < len(proto_init[0]):
###################################################################
updated_inputs = {"Input temperature [K]": spm_temperature}
vlims_ok = tmp_manager._step(step, proto_init[0], term_init[0], type_init[0], updated_inputs, skip_vcheck)
###################################################################
# Apply Heat Sources
Q_tot = tmp_manager.output[Qid, step, :]
Q = get_cc_power_loss(net, netlist)
# To do - Get cc heat from netlist
# Q_ohm_cc = net.interpolate_data("pore.cc_power_loss")[res_Ts]
# Q_ohm_cc /= net["throat.volume"][res_Ts]
# key = "Volume-averaged Ohmic heating CC [W.m-3]"
# vh[key][outer_step, :] = Q_ohm_cc[sorted_res_Ts]
Q[res_Ts] += Q_tot
jellybamm.apply_heat_source_lp(project, Q)
# Calculate Global Temperature
jellybamm.run_step_transient(
project, dim_time_step, T0, cp, rho, thermal_third
)
# Interpolate the node temperatures for the SPMs
spm_temperature = phase.interpolate_data("pore.temperature")[res_Ts]
# T_non_dim_spm = fT_non_dim(parameter_values, spm_temperature)
###################################################################
if vlims_ok:
step += 1
pbar.update(1)
temp_Ri = np.array(netlist.loc[tmp_manager.Ri_map].value)
# plt.scatter(np.arange(len(temp_Ri)), temp_Ri, label=str(step))
else:
break
# plt.legend()
tmp_manager.step = step
return netlist, project, phase, spm_temperature
toc = ticker.time()
lp.logger.notice("Initial step solve finished")
lp.logger.notice("Total stepping time " + str(np.around(toc - tic, 3)) + "s")
lp.logger.notice(
"Time per step " + str(np.around((toc - tic) / tmp_manager.Nsteps, 3)) + "s"
)


netlist, project, phase, spm_temperature = initialize_simulation(I_app, netlist, project, phase, spm_temperature)
###########################################################################
# Real Solve
###########################################################################
Expand All @@ -221,26 +228,35 @@ def run_simulation_lp(parameter_values, experiment, initial_soc, project):
Qid = np.argwhere(np.asarray(manager.variable_names) == Qvar).flatten()[0]
netlist["power_loss"] = 0.0
manager.global_step = 0
step_previous_Iapp = 0
print(types, terms)
for ps, step_protocol in enumerate(protos):
step_termination = terms[ps]
step_type = steps[ps]
step_type = types[ps]
if step_termination == []:
step_termination = 0.0
# normally would run this:


# do reinitialization if switching directions
if step_protocol[0] != step_previous_Iapp and step_protocol[0] != 0:
netlist, project, phase, spm_temperature = initialize_simulation(step_protocol[0], netlist, project, phase, spm_temperature)
manager.netlist = netlist
step_previous_Iapp = step_protocol[-1]

## Now Solve
tic = ticker.time()
lp.logger.notice("Starting step solve")
lp.logger.notice(f"Starting step solve for step {ps}")

vlims_ok = True
skip_vcheck = True
with tqdm(total=len(step_protocol), desc="Stepping simulation") as pbar:
with tqdm(total=len(step_protocol), desc=f"Stepping simulation ({ps+1}/{len(protos)})") as pbar:
step = 0
while step < len(step_protocol):
###################################################################
updated_inputs = {"Input temperature [K]": spm_temperature}
vlims_ok = manager._step(step, step_protocol, step_termination, step_type, updated_inputs, skip_vcheck)
skip_vcheck = False
if step > 5:
skip_vcheck = False
###################################################################
# Apply Heat Sources
Q_tot = manager.output[Qid, step, :]
Expand All @@ -264,6 +280,7 @@ def run_simulation_lp(parameter_values, experiment, initial_soc, project):
step += 1
pbar.update(1)
else:
# This break statement will only break out of the innermost while loop
break
manager.step = step
toc = ticker.time()
Expand All @@ -273,7 +290,7 @@ def run_simulation_lp(parameter_values, experiment, initial_soc, project):
"Time per step " + str(np.around((toc - tic) / manager.Nsteps, 3)) + "s"
)

print("*" * 30)
print("ECM Sim time", ticker.time() - st)
print("*" * 30)
return project, manager.step_output()
print("*" * 30)
print("Step Sim time", ticker.time() - st)
print("*" * 30)
return project, manager.step_output()
104 changes: 101 additions & 3 deletions tests/unit/test_funcs.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@
wrk = op.Workspace()


class funcsTest(unittest.TestCase):
class OneStepExperimentTest(unittest.TestCase):
def _teardown(self, fpath):
# Delete Data files
pass
Expand All @@ -21,10 +21,108 @@ def _jellybamm_general(self, project):
I_app = 0.35
dt = 5
Nsteps = 3
hours = dt * Nsteps
seconds = dt * Nsteps
experiment = pybamm.Experiment(
[
f"Discharge at {I_app} A for {hours} seconds",
f"Discharge at {I_app} A for {seconds} seconds",
],
period=f"{dt} seconds",
)

# Parameter set
param = pybamm.ParameterValues("Chen2020")
# JellyBaMM discretises the spiral using the electrode height for spiral length
# This parameter set has the longer length set to the Electrode width
# We want to swap this round
param["Electrode width [m]"] = 0.08
# Passing None as initial_soc will take values from Parameter set and apply
# uniformly everywhere
initial_soc = None

# Run simulation
project, output = jellybamm.run_simulation_lp(
parameter_values=param,
experiment=experiment,
initial_soc=initial_soc,
project=project,
)
assert output is not None

def test_jellybamm_spiral(self):
wrk.clear()
Nlayers = 2
dtheta = 10
spacing = 195e-6 # To do should come from params
inner_r = 10 * spacing
pos_tabs = [-1]
neg_tabs = [0]
length_3d = 0.08
tesla_tabs = False
# OpenPNM project
project, arc_edges = jellybamm.make_spiral_net(
Nlayers, dtheta, spacing, inner_r, pos_tabs, neg_tabs, length_3d, tesla_tabs
)
self._jellybamm_general(project)

def test_jellybamm_spiral_tesla(self):
wrk.clear()
Nlayers = 2
dtheta = 10
spacing = 195e-6 # To do should come from params
inner_r = 10 * spacing
pos_tabs = [-1]
neg_tabs = [0]
length_3d = 0.08
tesla_tabs = True
# OpenPNM project
project, arc_edges = jellybamm.make_spiral_net(
Nlayers, dtheta, spacing, inner_r, pos_tabs, neg_tabs, length_3d, tesla_tabs
)
self._jellybamm_general(project)

def test_jellybamm_tomo(self):
wrk.clear()
# Geometry of spiral
tomo_pnm = "spider_net.pnm"
dtheta = 10
spacing = 195e-6
length_3d = 0.08
pos_tabs = [-1]
neg_tabs = [0]
# OpenPNM project
project, arc_edges = jellybamm.make_tomo_net(
tomo_pnm, dtheta, spacing, length_3d, pos_tabs, neg_tabs
)
self._jellybamm_general(project)

def test_jellybamm_1d(self):
wrk.clear()
# Geometry of 1D mesh
Nunit = 100
spacing = 0.01
pos_tabs = [-1]
neg_tabs = [0]
# OpenPNM project
project, arc_edges = jellybamm.make_1D_net(Nunit, spacing, pos_tabs, neg_tabs)
self._jellybamm_general(project)


class MultiStepExperimentTest(unittest.TestCase):
def _teardown(self, fpath):
# Delete Data files
pass

def _jellybamm_general(self, project):
# Experiment
I_app = 0.35
dt = 5
Nsteps = 3
seconds = dt * Nsteps
experiment = pybamm.Experiment(
[
(f"Discharge at {I_app} A for {seconds} seconds"),
(f"Rest for {seconds} seconds"),
(f"Charge at {I_app} A for {seconds} seconds"),
],
period=f"{dt} seconds",
)
Expand Down

0 comments on commit 70a647e

Please sign in to comment.