Tidy up verification script and cell specification

cells/LGM50.hpp

@@ -21,7 +21,8 @@ const real_t kn_dim = 6.48e-7;      // Negative electrode reaction rate [A m^{2.
 const real_t Dn = 3.3e-14;          // Diffusion coefficient, negative electrode [m^2/s]
 const real_t An = 3 * eps_n_s / rn; // Negative electrode area [m^{-1}]
 
-const real_t sig_n = 215.; // Negative electrode conductivity [S/m]
+const real_t sig_n = 215.;  // Negative electrode conductivity [S/m]
+const real_t brugg_n = 1.5; // Negative electrode Bruggeman coefficient
 
 inline const real_t
 Un(real_t cs) // Negative electrode open circuit potential [V]
@@ -45,7 +46,8 @@ const real_t kp_dim = 3.42e-6;      // Positive electrode reaction rate [A m^{2.
 const real_t Dp = 4.0e-15;          // Diffusion coefficient, positive electrode [m^2/s]
 const real_t Ap = 3 * eps_p_s / rp; // Positive electrode area [m^2]
 
-const real_t sig_p = 0.18; // Positive electrode conductivity [S/m]
+const real_t sig_p = 0.18;  // Positive electrode conductivity [S/m]
+const real_t brugg_p = 1.5; // Positive electrode Bruggeman coefficient
 
 inline const real_t
 Up(real_t cs) // Positive electrode open circuit potential [V]
@@ -55,29 +57,29 @@ Up(real_t cs) // Positive electrode open circuit potential [V]
 }
 
 // Separator.
-const real_t ls = 12e-6;     // Separator thickness [m]
-const real_t tplus = 0.2594; // Transference number of Li ions
-const real_t eps_s = 0.47;   // Separator electrolyte volume fraction
+const real_t ls = 12e-6;    // Separator thickness [m]
+const real_t eps_s = 0.47;  // Separator electrolyte volume fraction
+const real_t brugg_s = 1.5; // Separator Bruggeman coefficient
 
 // Electrolyte.
-const real_t ce0 = 1000.0; // Initial electrolyte concentration [mol/(m^3)]
+const real_t tplus = 0.2594; // Transference number of Li ions
+const real_t ce0 = 1000.0;   // Initial electrolyte concentration [mol/(m^3)]
 inline const real_t
 De(real_t ce) // Diffusivity of Li ions in the electrolyte [m^2/s]
 {
   return 8.794e-11 * pow(ce / 1000, 2) - 3.972e-10 * (ce / 1000) + 4.862e-10;
 }
 inline const real_t
-kappa(real_t x) // Electronic conductivity [S/m]
+kappa(real_t ce) // Electronic conductivity [S/m]
 {
-  return 0.1297 * pow(x / 1000, 3) - 2.51 * pow(x / 1000, 1.5) + 3.329 * (x / 1000);
+  return 0.1297 * pow(ce / 1000, 3) - 2.51 * pow(ce / 1000, 1.5) + 3.329 * (ce / 1000);
 }
 
 // Cell parameters.
-const real_t llayer = 1.58;               // Wound layer length [m]
-const real_t wlayer = 6.5e-2;             // Wound layer width [m]
-const real_t cell_area = wlayer * llayer; // Wound layer area (one layer for wound cells) [m^2]
+const real_t llayer = 1.58;                        // Layer length [m]
+const real_t wlayer = 6.5e-2;                      // Layer width [m]
+const unsigned nlayer = 1;                         // Number of layers (one layer for wound cells)
+const real_t cell_area = nlayer * wlayer * llayer; // Total area [m^2]
 
 const real_t I_typ = 5.0; // Or I1C in Jubat. Reference total current [A]
-
-const real_t brugg = 1.5; // Currently the same in each region
 }

parameters/constants.hpp

@@ -73,9 +73,9 @@ const real_t De_scale = L * L / te;
 
 // Transport efficiency (inverse MacMullin number). This is B(x) in Planella, and is absorbed into
 // the definition of kappa_ne/kappa_pe/kappa_sp in JuBat.
-const real_t BPE = pow(eps_p, brugg);
-const real_t BNE = pow(eps_n, brugg);
-const real_t BSEP = pow(eps_s, brugg);
+const real_t BPE = pow(eps_p, brugg_p);
+const real_t BNE = pow(eps_n, brugg_n);
+const real_t BSEP = pow(eps_s, brugg_s);
 
 const real_t j_scale = I_typ / a0 / L / cell_area;
 
@@ -140,9 +140,9 @@ DE(real_t ce)
   return De(ce * ce_scale) / De_scale;
 }
 inline const real_t
-Kappa(real_t x)
+Kappa(real_t ce)
 {
-  return kappa(x * ce_scale) / kappa_scale;
+  return kappa(ce * ce_scale) / kappa_scale;
 }
 
 const real_t KS = Kappa(CE0); // / kappa_scale; // Scaled electrolyte conductivity.

parameters/settings.hpp

@@ -9,9 +9,8 @@ extern const bool P2D;  // Whether running P2D model
 extern const unsigned NNE;  // Number of elements in the Negative Electrode
 extern const unsigned NSEP; // Number of elements in the Separator
 extern const unsigned NPE;  // Number of elements in the Positive Electrode
-extern const unsigned
-    NX; // Number of elements in the X-dimension (i.e Electrolye) (Sum of the above three)
-extern const unsigned NR; // Number of elements in the R-dimension (i.e Particle)
+extern const unsigned NX;   // Number of elements in the X-dimension (i.e Electrolye) (sum of above)
+extern const unsigned NR;   // Number of elements in the R-dimension (i.e Particle)
 
 extern const unsigned NNEPAR; // Number of Negative Electrode PARticle
 extern const unsigned NPEPAR; // Number of Positive Electrode PARticle

validation/run_pybamm_and_compare_to_mfem.py

@@ -1,6 +1,4 @@
-import pandas as pd
 import matplotlib.pyplot as plt
-import numpy as np
 import os
 import pybamm
 import subprocess
@@ -49,7 +47,7 @@ def easyplot(x, y, colour, linestyle, fig, label, marker_indx):
 # Function to run Batree and extract results from output.
 def run_batree(mfem_executable,sim_type):
 
-    print(''.join(["Running ", sim_type, " simulation in batree..."]))
+    print(''.join(["Running ", sim_type, " simulation in Batree..."]))
     cmd = [mfem_executable, "-m", sim_type]
     result = subprocess.run(cmd, capture_output=True, text=True)
 
@@ -59,8 +57,8 @@ def run_batree(mfem_executable,sim_type):
     for line in result.stdout.splitlines():
         parts = line.split()
 
-        if len(parts) == 3 and parts[0].isdigit():
-            step, time, voltage = parts
+        if len(parts) == 4 and parts[0].isdigit():
+            _, time, voltage, _ = parts
             times.append(float(time))
             voltages.append(float(voltage))
 
@@ -84,17 +82,17 @@ def run_pybamm(model, parameter_values):
 if PLOT_SPM and PLOT_MFEM:
 
     time, voltage = run_batree(mfem_executable, "SPM")
-    easyplot(time, voltage, blue, "-o", 1, "SPM (MFEM)", mfem_indx)
+    easyplot(time, voltage, blue, "-o", 1, "SPM (Batree)", mfem_indx)
 
 if PLOT_SPMe and PLOT_MFEM:
 
     time, voltage = run_batree(mfem_executable, "SPMe")
-    easyplot(time, voltage, red, "-o", 1, "SPMe (MFEM)", mfem_indx)
+    easyplot(time, voltage, red, "-o", 1, "SPMe (Batree)", mfem_indx)
 
 if PLOT_DFN and PLOT_MFEM:
 
     time, voltage = run_batree(mfem_executable, "P2D")
-    easyplot(time, voltage, black, "-o", 1, "DFN (MFEM)", mfem_indx)
+    easyplot(time, voltage, black, "-o", 1, "DFN (Batree)", mfem_indx)
 
 if PLOT_SPM and PLOT_PYBAMM: