Updates for 24 1

.gitignore

@@ -132,3 +132,5 @@ dmypy.json
 
 # Data
 data/
+
+*.pdf

examples/compare_methods.py

@@ -38,7 +38,7 @@ def current_function(t):
     # Extract final two periods of the solution
     time = sol["Time [s]"].entries[-3 * samples_per_period - 1 :]
     current = sol["Current [A]"].entries[-3 * samples_per_period - 1 :]
-    voltage = sol["Terminal voltage [V]"].entries[-3 * samples_per_period - 1 :]
+    voltage = sol["Voltage [V]"].entries[-3 * samples_per_period - 1 :]
     # FFT
     current_fft = fft(current)
     voltage_fft = fft(voltage)
@@ -52,7 +52,7 @@ def current_function(t):
 print("Time domain method: ", time_elapsed, "s")
 
 # Frequency domain
-methods = ["direct", "prebicgstab"]
+methods = ["direct"]
 impedances_freqs = []
 for method in methods:
     start_time = timer.time()
@@ -72,5 +72,4 @@ def current_function(t):
     )
 ax.legend()
 plt.suptitle(f"{model.name}")
-plt.savefig(f"figures/{model.name}_time_vs_freq.pdf", dpi=300)
 plt.show()

examples/compare_models.py

@@ -4,7 +4,13 @@
 
 # Load models and parameters
 models = [
-    pybamm.lithium_ion.SPM(options={"surface form": "differential"}, name="SPM"),
+    pybamm.lithium_ion.SPM(
+        options={
+            "surface form": "differential",
+        },
+        name="SPM",
+    ),
+    pybamm.lithium_ion.SPMe(options={"surface form": "differential"}, name="SPMe"),
     pybamm.lithium_ion.DFN(options={"surface form": "differential"}, name="DFN"),
     pybamm.lithium_ion.SPM(
         {
@@ -55,5 +61,4 @@
         impedances[i], ax=ax, linestyle="-", label=f"{model.name}", alpha=0.7
     )
 ax.legend()
-plt.savefig("figures/compare_models.pdf", dpi=300)
 plt.show()

pbeis/eis_simulation.py

@@ -76,9 +76,11 @@ def __init__(
         # final entry by construction
         self.b = np.zeros_like(self.y0)
         self.b[-1] = -1
-        # Store time and current scales
-        self.timescale = self.built_model.timescale_eval
-        self.current_scale = sim.parameter_values.evaluate(model.param.I_typ)
+        # Get scale factor for the impedance (PyBaMM model may have set scales for the
+        # voltage and current variables)
+        V_scale = getattr(self.model.variables["Voltage [V]"], "scale", 1)
+        I_scale = getattr(self.model.variables["Current [A]"], "scale", 1)
+        self.z_scale = parameter_values.evaluate(V_scale / I_scale)
 
         # Set setup time
         self.set_up_time = timer.time()
@@ -102,25 +104,24 @@ def set_up_model_for_eis(self, model):
         new_model = model.new_copy()
 
         # Create a voltage variable
-        V_cell = pybamm.Variable("Terminal voltage variable")
-        new_model.variables["Terminal voltage variable"] = V_cell
-        V = new_model.variables["Terminal voltage [V]"]
+        V_cell = pybamm.Variable("Voltage variable [V]")
+        new_model.variables["Voltage variable [V]"] = V_cell
+        V = new_model.variables["Voltage [V]"]
         # Add an algebraic equation for the voltage variable
         new_model.algebraic[V_cell] = V_cell - V
-        new_model.initial_conditions[V_cell] = (
-            new_model.param.p.U_ref - new_model.param.n.U_ref
-        )
+        new_model.initial_conditions[V_cell] = new_model.param.ocv_init
 
         # Now make current density a variable
-        # To do so, we replace all instances of the current density in the
-        # model with a current density variable, which is obtained from the
-        # FunctionControl submodel
+        # To do so, we replace all instances of the current in the model with a current
+        # variable, which is obtained from the FunctionControl submodel
 
         # Create the FunctionControl submodel and extract variables
         external_circuit_variables = pybamm.external_circuit.FunctionControl(
             model.param, None, model.options, control="algebraic"
         ).get_fundamental_variables()
 
+        # TODO: remove SymbolReplacer and use PyBaMM's "set up for experiment"
+
         # Perform the replacement
         symbol_replacement_map = {
             new_model.variables[name]: variable
@@ -133,15 +134,15 @@ def set_up_model_for_eis(self, model):
         )
         replacer.process_model(new_model, inplace=True)
 
-        # Add an algebraic equation for the current density variable
+        # Add an algebraic equation for the current variable
         # External circuit submodels are always equations on the current
-        i_cell = new_model.variables["Current density variable"]
+        I_var = new_model.variables["Current variable [A]"]
         I = new_model.variables["Current [A]"]
         I_applied = pybamm.FunctionParameter(
-            "Current function [A]", {"Time [s]": pybamm.t * new_model.param.timescale}
+            "Current function [A]", {"Time [s]": pybamm.t}
         )
-        new_model.algebraic[i_cell] = I - I_applied
-        new_model.initial_conditions[i_cell] = 0
+        new_model.algebraic[I_var] = I - I_applied
+        new_model.initial_conditions[I_var] = 0
 
         pybamm.logger.info("Finish setting up {} for EIS".format(self.model_name))
 
@@ -178,7 +179,7 @@ def solve(self, frequencies, method="direct"):
         if method == "direct":
             zs = []
             for frequency in frequencies:
-                A = 1.0j * 2 * np.pi * frequency * self.timescale * self.M - self.J
+                A = 1.0j * 2 * np.pi * frequency * self.M - self.J
                 lu = splu(A)
                 x = lu.solve(self.b)
                 # The model is set up such that the voltage is the penultimate
@@ -193,9 +194,7 @@ def solve(self, frequencies, method="direct"):
                 f"but is '{method}'",
             )
 
-        # Note: the current density variable is dimensionless so we need
-        # to scale by the current scale from the model to get true impedance
-        self.solution = np.array(zs) / self.current_scale
+        self.solution = np.array(zs) * self.z_scale
 
         # Store solve time as an attribute
         self.solve_time = timer.time()
@@ -246,7 +245,7 @@ def iterative_method(self, frequencies, method="prebicgstab"):
             num_iters = 0
 
             # Construct the matrix A(frequency)
-            A = 1.0j * 2 * np.pi * frequency * self.timescale * self.M - self.J
+            A = 1.0j * 2 * np.pi * frequency * self.M - self.J
 
             def callback(xk):
                 """
@@ -286,7 +285,7 @@ def callback(xk):
             z = -sol[-2][0] / sol[-1][0]
             zs.append(z)
 
-        return zs
+        return self.zs
 
     def nyquist_plot(self, ax=None, marker="o", linestyle="None", **kwargs):
         """

pbeis/plotting.py

@@ -36,9 +36,13 @@ def nyquist_plot(data, ax=None, marker="o", linestyle="None", **kwargs):
         show = False
 
     ax.plot(data.real, -data.imag, marker=marker, linestyle=linestyle, **kwargs)
+    _, xmax = ax.get_xlim()
+    _, ymax = ax.get_ylim()
+    axmax = max(xmax, ymax)
+    plt.axis([0, axmax, 0, axmax])
+    plt.gca().set_aspect("equal", adjustable="box")
     ax.set_xlabel(r"$Z_\mathrm{Re}$ [Ohm]")
     ax.set_ylabel(r"$-Z_\mathrm{Im}$ [Ohm]")
-
     if show:
         plt.show()
 

pbeis/utils.py

@@ -74,9 +74,9 @@ def process_model(self, unprocessed_model, inplace=True):
                 "Replacing symbols in {!r} (initial conditions)".format(variable)
             )
             if self.process_initial_conditions:
-                new_initial_conditions[
-                    self.process_symbol(variable)
-                ] = self.process_symbol(equation)
+                new_initial_conditions[self.process_symbol(variable)] = (
+                    self.process_symbol(equation)
+                )
             else:
                 new_initial_conditions[self.process_symbol(variable)] = equation
         model.initial_conditions = new_initial_conditions
@@ -101,20 +101,6 @@ def process_model(self, unprocessed_model, inplace=True):
             )
         model.events = new_events
 
-        # Set external variables
-        model.external_variables = [
-            self.process_symbol(var) for var in unprocessed_model.external_variables
-        ]
-
-        # Process timescale
-        model._timescale = self.process_symbol(unprocessed_model.timescale)
-
-        # Process length scales
-        new_length_scales = {}
-        for domain, scale in unprocessed_model.length_scales.items():
-            new_length_scales[domain] = self.process_symbol(scale)
-        model._length_scales = new_length_scales
-
         pybamm.logger.info("Finish replacing symbols in {}".format(model.name))
 
         return model

setup.py

@@ -15,5 +15,8 @@
     author="Rishit Dhoot & Robert Timms",
     author_email="timms@maths.ox.ac.uk",
     license="LICENSE",
-    install_requires=["pybamm == 22.10", "matplotlib"],
+    install_requires=[
+        "PyBaMM @ git+https://github.com/pybamm-team/PyBaMM.git@c29e736",
+        "matplotlib",
+    ],
 )