Skip to content

Commit

Permalink
#1129 running more julia tests
Browse files Browse the repository at this point in the history
  • Loading branch information
valentinsulzer committed Nov 11, 2020
1 parent bec5ce3 commit f3e4936
Show file tree
Hide file tree
Showing 2 changed files with 44 additions and 37 deletions.
41 changes: 21 additions & 20 deletions examples/scripts/DFN.py
Original file line number Diff line number Diff line change
Expand Up @@ -20,9 +20,10 @@
param.process_geometry(geometry)

# set mesh
# var = pybamm.standard_spatial_vars
# var_pts = {var.x_n: 30, var.x_s: 30, var.x_p: 30, var.r_n: 10, var.r_p: 10}
mesh = pybamm.Mesh(geometry, model.default_submesh_types, model.default_var_pts)
var = pybamm.standard_spatial_vars
var_pts = {var.x_n: 30, var.x_s: 30, var.x_p: 30, var.r_n: 10, var.r_p: 10}
# var_pts = model.default_var_pts
mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts)

# discretise model
disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
Expand All @@ -33,22 +34,22 @@
solver = pybamm.CasadiSolver(mode="fast", atol=1e-6, rtol=1e-6)
solution = solver.solve(model, t_eval)
solution = solver.solve(model, t_eval)
print(solution.integration_time)
print(pybamm.Timer().format(solution.integration_time))

# plot
plot = pybamm.QuickPlot(
solution,
[
"Negative particle concentration [mol.m-3]",
"Electrolyte concentration [mol.m-3]",
"Positive particle concentration [mol.m-3]",
"Current [A]",
"Negative electrode potential [V]",
"Electrolyte potential [V]",
"Positive electrode potential [V]",
"Terminal voltage [V]",
],
time_unit="seconds",
spatial_unit="um",
)
plot.dynamic_plot()
# plot = pybamm.QuickPlot(
# solution,
# [
# "Negative particle concentration [mol.m-3]",
# "Electrolyte concentration [mol.m-3]",
# "Positive particle concentration [mol.m-3]",
# "Current [A]",
# "Negative electrode potential [V]",
# "Electrolyte potential [V]",
# "Positive electrode potential [V]",
# "Terminal voltage [V]",
# ],
# time_unit="seconds",
# spatial_unit="um",
# )
# plot.dynamic_plot()
40 changes: 23 additions & 17 deletions tests/unit/test_expression_tree/test_operations/quick_julia_test.py
Original file line number Diff line number Diff line change
Expand Up @@ -17,33 +17,37 @@
from julia import Main


model = pybamm.lithium_ion.DFN()
sim = pybamm.Simulation(model, solver=pybamm.CasadiSolver(mode="fast"))
model = pybamm.lithium_ion.SPM()
# var = pybamm.standard_spatial_vars
# var_pts = {var.x_n: 30, var.x_s: 30, var.x_p: 30, var.r_n: 10, var.r_p: 10}
var_pts = model.default_var_pts
sim = pybamm.Simulation(model, solver=pybamm.CasadiSolver(mode="fast"), var_pts=var_pts)
sim.solve([0, 3600])
param = model.default_parameter_values
timescale = param.evaluate(model.timescale)
sol = sim.solve(np.linspace(0, 0.15 * timescale, 100))
print(sol.y[:, -1])
print(sol.integration_time)
expr = pybamm.NumpyConcatenation(
sim.built_model.concatenated_rhs, sim.built_model.concatenated_algebraic
).simplify()
# expr = sim.built_model.concatenated_rhs.simplify()
# expr = pybamm.NumpyConcatenation(
# sim.built_model.concatenated_rhs, sim.built_model.concatenated_algebraic
# ).simplify()
expr = sim.built_model.concatenated_rhs.simplify()

evaluator_str = pybamm.get_julia_function(expr)
n_rhs = sim.built_model.concatenated_rhs.size
n_alg = sim.built_model.concatenated_algebraic.size
# np.set_printoptions(
# threshold=max(
# np.get_printoptions()["threshold"],
# n_rhs + n_alg,
# )
# threshold=100
# # max(
# # np.get_printoptions()["threshold"],
# # n_rhs + n_alg,
# # )
# )
with open("tmp.txt", "w") as f:
with open("tmp_SPM_30.txt", "w") as f:
f.write(evaluator_str + "\n\n")
f.write(f"u0 = {np.array2string(sol.model.y0, separator=',')}\n")
f.write(f"du0 = zeros({n_rhs + n_alg})\n")
f.write(f"differential_vars=[ones({n_rhs});zeros({n_alg})]\n")
f.write(f"u0 = {np.array2string(sol.model.y0.full().flatten(), separator=',')}\n")
# f.write(f"du0 = zeros({n_rhs + n_alg})\n")
# f.write(f"differential_vars=[ones({n_rhs});zeros({n_alg})]\n")

# expr2 = sim.built_model.variables["Terminal voltage [V]"]
# evaluator_str2 = pybamm.get_julia_function(expr2)
Expand All @@ -52,10 +56,12 @@

Main.eval(evaluator_str)
Main.dy = np.zeros(n_rhs + n_alg)
Main.y = sol.model.y0
Main.y = sol.model.y0.full().flatten()
Main.eval("f(dy,y,0,0)")
print(Main.dy)
expr.evaluate(y=sol.model.y0)
# print(Main.dy)
# expr.evaluate(y=sol.model.y0)

print(Main.dy - expr.evaluate(y=Main.y).T)
# print(expr.evaluate(y=Main.y))
# # test something with a heaviside
# a = pybamm.Vector([1, 2])
Expand Down

0 comments on commit f3e4936

Please sign in to comment.