pybamm-team#744-added Functions for particle size distribution

.gitignore

@@ -4,6 +4,7 @@
 *.png
 /local/
 *.DS_Store
+*.mat
 
 # don't ignore important .txt files
 !requirements*

examples/scripts/compare_lithium_ion_particle_distribution.py

@@ -0,0 +1,64 @@
+#
+# Compare lithium-ion battery models
+#
+import argparse
+import numpy as np
+import pybamm
+
+parser = argparse.ArgumentParser()
+parser.add_argument(
+    "--debug", action="store_true", help="Set logging level to 'DEBUG'."
+)
+args = parser.parse_args()
+if args.debug:
+    pybamm.set_logging_level("DEBUG")
+else:
+    pybamm.set_logging_level("INFO")
+
+# load models
+options = {"thermal": "isothermal"}
+models = [
+    pybamm.lithium_ion.DFN(options,name = 'standard DFN'),
+    pybamm.lithium_ion.DFN(options,name = 'particle DFN'),
+]
+
+
+# load parameter values and process models and geometry
+params = [models[0].default_parameter_values, models[1].default_parameter_values]
+params[0]["Typical current [A]"] = 1.0
+params[0].process_model(models[0])
+
+
+params[1]["Typical current [A]"] = 1.0
+
+def negative_distribution(x):
+    return 1 + x
+def positive_distribution(x):
+    return 1 + (x-(1-models[1].param.l_p))
+
+params[1]["Negative particle distribution"] = negative_distribution
+params[1]["Positive particle distribution"] = positive_distribution
+params[1].process_model(models[1])
+
+# set mesh
+var = pybamm.standard_spatial_vars
+var_pts = {var.x_n: 10, var.x_s: 10, var.x_p: 10, var.r_n: 5, var.r_p: 5}
+
+# discretise models
+for param, model in zip(params, models):
+    # create geometry
+    geometry = model.default_geometry
+    param.process_geometry(geometry)
+    mesh = pybamm.Mesh(geometry, models[-1].default_submesh_types, var_pts)
+    disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
+    disc.process_model(model)
+
+# solve model
+solutions = [None] * len(models)
+t_eval = np.linspace(0, 0.3, 100)
+for i, model in enumerate(models):
+    solutions[i] = model.default_solver.solve(model, t_eval)
+
+# plot
+plot = pybamm.QuickPlot(models, mesh, solutions)
+plot.dynamic_plot()

input/parameters/lithium-ion/anodes/graphite_mcmb2528_Marquis2019/negative_particle_distribution.py

@@ -0,0 +1,6 @@
+import pybamm 
+import numpy as np
+
+def negative_particle_distribution(x):
+
+    return pybamm.Function(np.ones_like,x)

input/parameters/lithium-ion/anodes/graphite_mcmb2528_Marquis2019/parameters.csv

@@ -11,6 +11,7 @@ Negative electrode OCP [V],[function]graphite_mcmb2528_ocp_Dualfoil1998,
 Negative electrode porosity,0.3,Scott Moura FastDFN,electrolyte volume fraction
 Negative electrode active material volume fraction,0.7,,assuming zero binder volume fraction
 Negative particle radius [m],1E-05,Scott Moura FastDFN,
+Negative particle distribution,[function]negative_particle_distribution,,
 Negative electrode surface area density [m-1],180000,Scott Moura FastDFN,
 Negative electrode Bruggeman coefficient (electrolyte),1.5,Scott Moura FastDFN,
 Negative electrode Bruggeman coefficient (electrode),1.5,Scott Moura FastDFN,

input/parameters/lithium-ion/cathodes/lico2_Marquis2019/parameters.csv

@@ -11,6 +11,7 @@ Positive electrode OCP [V],[function]lico2_ocp_Dualfoil1998,
 Positive electrode porosity,0.3,Scott Moura FastDFN,electrolyte volume fraction
 Positive electrode active material volume fraction,0.7,,assuming zero binder volume fraction
 Positive particle radius [m],1E-05,Scott Moura FastDFN,
+Positive particle distribution,[function]positive_particle_distribution,,
 Positive electrode surface area density [m-1],150000,Scott Moura FastDFN,
 Positive electrode Bruggeman coefficient (electrolyte),1.5,Scott Moura FastDFN,
 Positive electrode Bruggeman coefficient (electrode),1.5,Scott Moura FastDFN,

input/parameters/lithium-ion/cathodes/lico2_Marquis2019/positive_particle_distribution.py

@@ -0,0 +1,6 @@
+import pybamm 
+import numpy as np
+
+def positive_particle_distribution(x):
+
+    return pybamm.Function(np.ones_like,x)

pybamm/expression_tree/operations/convert_to_casadi.py

@@ -75,6 +75,9 @@ def _convert(self, symbol, t, y):
                 return casadi.mmax(*converted_children)
             elif symbol.function == np.abs:
                 return casadi.fabs(*converted_children)
+            elif symbol.function == np.ones_like:
+                len_child = converted_children[0].shape[0]
+                return casadi.MX.ones(len_child)
             elif isinstance(symbol.function, (PchipInterpolator, CubicSpline)):
                 return casadi.interpolant("LUT", "bspline", [symbol.x], symbol.y)(
                     *converted_children

pybamm/models/submodels/particle/fickian/fickian_many_particles.py

@@ -55,12 +55,12 @@ def set_rhs(self, variables):
 
         if self.domain == "Negative":
             x = pybamm.standard_spatial_vars.x_n
-            R = pybamm.PrimaryBroadcast(1 + x / 10, 'negative particle',)
+            R = pybamm.PrimaryBroadcast(pybamm.FunctionParameter('Negative particle distribution',x), 'negative particle',)
             self.rhs = {c: -(1 / (R**2 * self.param.C_n)) * pybamm.div(N)}
 
         elif self.domain == "Positive":
             x = pybamm.standard_spatial_vars.x_p
-            R = pybamm.PrimaryBroadcast(1 + x / 10, 'positive particle',)
+            R = pybamm.PrimaryBroadcast(pybamm.FunctionParameter('Positive particle distribution', x), 'positive particle',)
             self.rhs = {c: -(1 / (R**2 * self.param.C_p)) * pybamm.div(N)}