tidy up SOC and make example more interesting + change solve terminat…

…ion to warning

examples/scripts/SPMe_SOC.py

@@ -1,110 +1,131 @@
+#
+#  Example to show the state of charge of a battery using the SPMe model
+#  Initial conditions are specified to start each electrode in 1/2 charged
+#  state. A charge and discharge are performed with current chosen to be
+#  1C rate when electrode dimensions are euqal.
+#  Coulomb counting is performed to calculate the capacity of the
+#  battery within the operating voltage limits and maximum particle concs.
+#  The anode thickenss is varied to highlight the importance of electrode
+#  sizing to enable full lithium utilization
+#
 import pybamm
 import numpy as np
 import matplotlib.pyplot as plt
 
 plt.close("all")
-pybamm.set_logging_level("INFO")
+pybamm.set_logging_level(30)
 
 factor = 6.38
-
-# Dimensions
-h = 0.137
-w = 0.207 / factor
-A = h * w
-l_n = 1e-4
+capacities = []
+specific_capacities = []
 l_p = 1e-4
-l_s = 2.5e-5
-l1d = (l_n + l_p + l_s)
-vol = h * w * l1d
-vol_cm3 = vol * 1e6
-
-tot_cap = 0.0
-tot_time = 0.0
-fig, axes = plt.subplots(1, 2, sharey=True)
-I_mag = 1.01 / factor
-for enum, I_app in enumerate([-1.0, 1.0]):
-    I_app *= I_mag
-    # load model
-    model = pybamm.lithium_ion.SPMe()
-    # create geometry
-    geometry = model.default_geometry
-    # load parameter values and process model and geometry
-    param = model.default_parameter_values
-
-    param.update(
-        {"Electrode height [m]": h,
-         "Electrode width [m]": w,
-         "Negative electrode thickness [m]": l_n,
-         "Positive electrode thickness [m]": l_p,
-         "Separator thickness [m]": l_s,
-         "Lower voltage cut-off [V]": 3.105,
-         "Upper voltage cut-off [V]": 4.7,
-         "Maximum concentration in negative electrode [mol.m-3]": 25000,
-         "Maximum concentration in positive electrode [mol.m-3]": 50000,
-         "Initial concentration in negative electrode [mol.m-3]": 12500,
-         "Initial concentration in positive electrode [mol.m-3]": 25000,
-         "Negative electrode surface area density [m-1]": 180000.0,
-         "Positive electrode surface area density [m-1]": 150000.0,
-         "Typical current [A]": I_app,
-         }
-    )
-
-    param.process_model(model)
-    param.process_geometry(geometry)
-    s_var = pybamm.standard_spatial_vars
-    var_pts = {s_var.x_n: 5, s_var.x_s: 5, s_var.x_p: 5,
-               s_var.r_n: 5, s_var.r_p: 10}
-    # set mesh
-    mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts)
-    # discretise model
-    disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
-    disc.process_model(model)
-    # solve model
-    t_eval = np.linspace(0, 0.2, 100)
-    sol = model.default_solver.solve(model, t_eval)
-
-    var = "Positive electrode average extent of lithiation"
-    xpext = pybamm.ProcessedVariable(model.variables[var],
-                                     sol.t, sol.y, mesh=mesh)
-    var = "Negative electrode average extent of lithiation"
-    xnext = pybamm.ProcessedVariable(model.variables[var],
-                                     sol.t, sol.y, mesh=mesh)
-    var = "X-averaged positive particle surface concentration"
-    xpsurf = pybamm.ProcessedVariable(model.variables[var],
-                                      sol.t, sol.y, mesh=mesh)
-    var = "X-averaged negative particle surface concentration"
-    xnsurf = pybamm.ProcessedVariable(model.variables[var],
-                                      sol.t, sol.y, mesh=mesh)
-    time = pybamm.ProcessedVariable(model.variables["Time [h]"],
-                                    sol.t, sol.y, mesh=mesh)
-
-    # Coulomb counting
-    time_hours = time(sol.t)
-    dc_time = np.around(time_hours[-1], 3)
-    # Capacity mAh
-    cap = np.absolute(I_app * 1000 * dc_time)
-    cap_time = np.absolute(I_app * 1000 * time_hours)
+thicknesses = np.linspace(1.0, 2.5, 11) * l_p
+for l_n in thicknesses:
+    e_ratio = np.around(l_n / l_p, 3)
+    # Dimensions
+    h = 0.137
+    w = 0.207 / factor
+    A = h * w
+    l_s = 2.5e-5
+    l1d = (l_n + l_p + l_s)
+    vol = h * w * l1d
+    vol_cm3 = vol * 1e6
+    tot_cap = 0.0
+    tot_time = 0.0
+    fig, axes = plt.subplots(1, 2, sharey=True)
+    I_mag = 1.01 / factor
+    print('*' * 30)
+    for enum, I_app in enumerate([-1.0, 1.0]):
+        I_app *= I_mag
+        # load model
+        model = pybamm.lithium_ion.SPMe()
+        # create geometry
+        geometry = model.default_geometry
+        # load parameter values and process model and geometry
+        param = model.default_parameter_values
+        param.update(
+            {"Electrode height [m]": h,
+             "Electrode width [m]": w,
+             "Negative electrode thickness [m]": l_n,
+             "Positive electrode thickness [m]": l_p,
+             "Separator thickness [m]": l_s,
+             "Lower voltage cut-off [V]": 2.8,
+             "Upper voltage cut-off [V]": 4.7,
+             "Maximum concentration in negative electrode [mol.m-3]": 25000,
+             "Maximum concentration in positive electrode [mol.m-3]": 50000,
+             "Initial concentration in negative electrode [mol.m-3]": 12500,
+             "Initial concentration in positive electrode [mol.m-3]": 25000,
+             "Negative electrode surface area density [m-1]": 180000.0,
+             "Positive electrode surface area density [m-1]": 150000.0,
+             "Typical current [A]": I_app,
+             }
+        )
+        param.process_model(model)
+        param.process_geometry(geometry)
+        s_var = pybamm.standard_spatial_vars
+        var_pts = {s_var.x_n: 5, s_var.x_s: 5, s_var.x_p: 5,
+                   s_var.r_n: 5, s_var.r_p: 10}
+        # set mesh
+        mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts)
+        # discretise model
+        disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
+        disc.process_model(model)
+        # solve model
+        t_eval = np.linspace(0, 0.2, 100)
+        sol = model.default_solver.solve(model, t_eval)
+        var = "Positive electrode average extent of lithiation"
+        xpext = pybamm.ProcessedVariable(model.variables[var],
+                                         sol.t, sol.y, mesh=mesh)
+        var = "Negative electrode average extent of lithiation"
+        xnext = pybamm.ProcessedVariable(model.variables[var],
+                                         sol.t, sol.y, mesh=mesh)
+        var = "X-averaged positive particle surface concentration"
+        xpsurf = pybamm.ProcessedVariable(model.variables[var],
+                                          sol.t, sol.y, mesh=mesh)
+        var = "X-averaged negative particle surface concentration"
+        xnsurf = pybamm.ProcessedVariable(model.variables[var],
+                                          sol.t, sol.y, mesh=mesh)
+        time = pybamm.ProcessedVariable(model.variables["Time [h]"],
+                                        sol.t, sol.y, mesh=mesh)
+        # Coulomb counting
+        time_hours = time(sol.t)
+        dc_time = np.around(time_hours[-1], 3)
+        # Capacity mAh
+        cap = np.absolute(I_app * 1000 * dc_time)
+        cap_time = np.absolute(I_app * 1000 * time_hours)
+        axes[enum].plot(cap_time,
+                        xnext(sol.t), 'r-', label='Average Neg')
+        axes[enum].plot(cap_time,
+                        xpext(sol.t), 'b-', label='Average Pos')
+        axes[enum].plot(cap_time,
+                        xnsurf(sol.t), 'r--', label='Surface Neg')
+        axes[enum].plot(cap_time,
+                        xpsurf(sol.t), 'b--', label='Surface Pos')
+        axes[enum].set_xlabel('Capacity [mAh]')
+        handles, labels = axes[enum].get_legend_handles_labels()
+        axes[enum].legend(handles, labels)
+        if I_app < 0.0:
+            axes[enum].set_ylabel('Extent of Lithiation, Elecrode Ratio: '
+                                  + str(e_ratio))
+            axes[enum].title.set_text('Charge')
+        else:
+            axes[enum].title.set_text('Discharge')
+        print('Applied Current', I_app, 'A', 'Time',
+              dc_time, 'hrs', 'Capacity', cap, 'mAh')
+        tot_cap += cap
+        tot_time += dc_time
 
-    axes[enum].plot(cap_time,
-                    xnext(sol.t), 'r-', label='Average Neg')
-    axes[enum].plot(cap_time,
-                    xpext(sol.t), 'b-', label='Average Pos')
-    axes[enum].plot(cap_time,
-                    xnsurf(sol.t), 'r--', label='Surface Neg')
-    axes[enum].plot(cap_time,
-                    xpsurf(sol.t), 'b--', label='Surface Pos')
-    axes[enum].set_xlabel('Capacity [mAh]')
-    plt.legend()
-    if I_app < 0.0:
-        axes[enum].set_ylabel('Extent of Lithiation')
-        axes[enum].title.set_text('Charge')
-    else:
-        axes[enum].title.set_text('Discharge')
-    print('Applied Current', I_app, 'A', 'Time',
-          dc_time, 'hrs', 'Capacity', cap, 'mAh')
-    tot_cap += cap
-    tot_time += dc_time
+    print('Anode : Cathode thickness', e_ratio)
+    print('Total Charge/Discharge Time', tot_time, 'hrs')
+    print('Total Capacity', np.around(tot_cap, 3), 'mAh')
+    specific_cap = np.around(tot_cap, 3) / vol_cm3
+    print('Total Capacity', specific_cap, 'mAh.cm-3')
+    capacities.append(tot_cap)
+    specific_capacities.append(specific_cap)
 
-print('Total Charge/Discharge Time', tot_time, 'hrs')
-print('Total Capacity', np.around(tot_cap, 3), 'mAh')
-print('Total Capacity', np.around(tot_cap, 3) / vol_cm3, 'mAh.cm-3')
+fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True)
+ax1.plot(thicknesses / l_p, capacities)
+ax2.plot(thicknesses / l_p, specific_capacities)
+ax1.set_ylabel('Capacity [mAh]')
+ax2.set_ylabel('Specific Capacity [mAh.cm-3]')
+ax2.set_xlabel('Anode : Cathode thickness')

pybamm/models/submodels/particle/base_particle.py

@@ -35,7 +35,7 @@ def _get_standard_concentration_variables(self, c_s, c_s_xav):
             c_scale = self.param.c_p_max
             active_volume = geo_param.a_p_dim * geo_param.R_p / 3
         c_s_r_av = pybamm.r_average(c_s_xav)
-        c_s_r_av_vol = active_volume * c_s_r_av * c_scale
+        c_s_r_av_vol = active_volume * c_s_r_av
         variables = {
             self.domain + " particle concentration": c_s,
             self.domain + " particle concentration [mol.m-3]": c_s * c_scale,
@@ -54,7 +54,9 @@ def _get_standard_concentration_variables(self, c_s, c_s_xav):
             + " particle surface concentration [mol.m-3]": c_scale * c_s_surf_av,
             self.domain + " electrode active volume fraction": active_volume,
             self.domain
-            + " electrode volume-averaged concentration [mol.m-3]": c_s_r_av_vol,
+            + " electrode volume-averaged concentration": c_s_r_av_vol,
+            self.domain + " electrode "
+            + "volume-averaged concentration [mol.m-3]": c_s_r_av_vol * c_scale,
             self.domain + " electrode average extent of lithiation": c_s_r_av
         }
 

pybamm/solvers/base_solver.py

@@ -64,7 +64,7 @@ def solve(self, model, t_eval):
         solution.total_time = timer.time() - start_time
         solution.set_up_time = set_up_time
 
-        pybamm.logger.info("Finish solving {} ({})".format(model.name, termination))
+        pybamm.logger.warning("Finish solving {} ({})".format(model.name, termination))
         pybamm.logger.info(
             "Set-up time: {}, Solve time: {}, Total time: {}".format(
                 timer.format(solution.set_up_time),