#617 exponential 1D mesh

examples/scripts/SPM_compare_particle_grid.py

@@ -1,75 +1,94 @@
-import pybamm
+#
+# Compare different discretisations in the particle
+#
+import argparse
 import numpy as np
+import pybamm
 import matplotlib.pyplot as plt
 
-pybamm.set_logging_level("INFO")
-
-# load model
-model = pybamm.lithium_ion.SPM()
+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")
 
-# create geometry
-geometry = model.default_geometry
+# load models
+models = [
+    pybamm.lithium_ion.SPM(name="Uniform mesh"),
+    pybamm.lithium_ion.SPM(name="Chebyshev mesh"),
+    pybamm.lithium_ion.SPM(name="Exponential mesh"),
+]
 
-# load parameter values and process model and geometry
-param = model.default_parameter_values
-param.process_model(model)
-param.process_geometry(geometry)
+# load parameter values and process models and geometry
+param = models[0].default_parameter_values
+param["Typical current [A]"] = 1.0
+for model in models:
+    param.process_model(model)
 
 # set mesh
-submesh_types = {
-    "negative electrode": pybamm.Uniform1DSubMesh,
-    "separator": pybamm.Uniform1DSubMesh,
-    "positive electrode": pybamm.Uniform1DSubMesh,
-    "negative particle": pybamm.Chebyshev1DSubMesh,
-    "positive particle": pybamm.Chebyshev1DSubMesh,
-    "current collector": pybamm.SubMesh0D,
-}
-var = pybamm.standard_spatial_vars
-var_pts = {var.x_n: 10, var.x_s: 10, var.x_p: 10, var.r_n: 9, var.r_p: 9}
-mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
-
-# discretise model
-disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
-disc.process_model(model)
+submesh_types = models[0].default_submesh_types
+particle_meshes = [
+    pybamm.Uniform1DSubMesh,
+    pybamm.Chebyshev1DSubMesh,
+    pybamm.RightExponential1DSubMesh,
+]
+meshes = [None] * len(models)
+# discretise models
+for i, model in enumerate(models):
+    # create geometry
+    geometry = model.default_geometry
+    param.process_geometry(geometry)
+    submesh_types["negative particle"] = particle_meshes[i]
+    submesh_types["positive particle"] = particle_meshes[i]
+    meshes[i] = pybamm.Mesh(geometry, submesh_types, model.default_var_pts)
+    disc = pybamm.Discretisation(meshes[i], model.default_spatial_methods)
+    disc.process_model(model)
 
 # solve model
-t_eval = np.linspace(0, 0.2, 100)
-solution = model.default_solver.solve(model, t_eval)
+solutions = [None] * len(models)
+t_eval = np.linspace(0, 0.17, 100)
+for i, model in enumerate(models):
+    solutions[i] = model.default_solver.solve(model, t_eval)
+
+# process particle concentration variables
+processed_variables = [None] * len(models)
+for i, solution in enumerate(solutions):
+    c_n = pybamm.ProcessedVariable(
+        models[i].variables["X-average negative particle concentration [mol.m-3]"],
+        solution.t,
+        solution.y,
+        mesh=meshes[i],
+    )
+    c_p = pybamm.ProcessedVariable(
+        models[i].variables["X-average positive particle concentration [mol.m-3]"],
+        solution.t,
+        solution.y,
+        mesh=meshes[i],
+    )
+    processed_variables[i] = {"c_n": c_n, "c_p": c_p}
+
 
 # plot
-r_n = mesh["negative particle"][0].edges
-c_n = pybamm.ProcessedVariable(
-    model.variables["X-average negative particle concentration [mol.m-3]"],
-    solution.t,
-    solution.y,
-    mesh=mesh,
-)
-r_p = mesh["positive particle"][0].edges
-c_p = pybamm.ProcessedVariable(
-    model.variables["X-average positive particle concentration [mol.m-3]"],
-    solution.t,
-    solution.y,
-    mesh=mesh,
-)
-import ipdb; ipdb.set_trace()
-fig, ax = plt.subplots(figsize=(15, 8))
-plt.tight_layout()
-plt.subplot(121)
-plt.plot(
-    r_n,
-    np.zeros_like(r_n),
-    "ro",
-    mesh["negative particle"][0].nodes,
-    c_n(t=0.1, r=mesh["negative particle"][0].nodes),
-    "b-",
-)
-plt.subplot(122)
-plt.plot(
-    r_p,
-    np.zeros_like(r_p),
-    "ro",
-    mesh["positive particle"][0].nodes,
-    c_p(t=0.1, r=mesh["positive particle"][0].nodes),
-    "b-",
-)
-plt.show()
+def plot(t):
+    fig, ax = plt.subplots(figsize=(15, 8))
+    plt.subplot(121)
+    plt.xlabel(r"$r_n$")
+    plt.ylabel("Negative particle concentration [mol.m-3]")
+    for i, vars in enumerate(processed_variables):
+        r_n = meshes[i]["negative particle"][0].nodes
+        neg_plot, = plt.plot(r_n, vars["c_n"](r=r_n, t=t), '-o', label=models[i].name)
+    plt.subplot(122)
+    plt.xlabel(r"$r_p$")
+    plt.ylabel("Positive particle concentration [mol.m-3]")
+    for i, vars in enumerate(processed_variables):
+        r_p = meshes[i]["positive particle"][0].nodes
+        pos_plot, = plt.plot(r_p, vars["c_p"](r=r_p, t=t), '-o', label=models[i].name)
+    plt.legend()
+    plt.show()
+
+
+plot(0.1)

pybamm/__init__.py

@@ -213,7 +213,14 @@ def version(formatted=False):
 from .discretisations.discretisation import Discretisation
 from .meshes.meshes import Mesh
 from .meshes.zero_dimensional_submesh import SubMesh0D
-from .meshes.one_dimensional_submeshes import SubMesh1D, Uniform1DSubMesh, Chebyshev1DSubMesh
+from .meshes.one_dimensional_submeshes import (
+    SubMesh1D,
+    Uniform1DSubMesh,
+    Chebyshev1DSubMesh,
+    SymmetricExponential1DSubMesh,
+    LeftExponential1DSubMesh,
+    RightExponential1DSubMesh,
+)
 from .meshes.scikit_fem_submeshes import Scikit2DSubMesh, have_scikit_fem
 
 #

pybamm/meshes/one_dimensional_submeshes.py

@@ -66,7 +66,6 @@ class Chebyshev1DSubMesh(SubMesh1D):
     .. math ::
      a < x_{1} < ... < x_{N} < b.
 
-
     Parameters
     ----------
     lims : dict
@@ -88,7 +87,7 @@ def __init__(self, lims, npts):
         npts = npts[spatial_var.id]
 
         # Create N Chebyshev nodes in the interval (a,b)
-        N = npts - 2
+        N = npts - 1
         ii = np.array(range(1, N + 1))
         a = spatial_lims["min"]
         b = spatial_lims["max"]
@@ -102,18 +101,19 @@ def __init__(self, lims, npts):
         super().__init__(edges, coord_sys=coord_sys)
 
 
-class ExponentialEdges1DSubMesh(SubMesh1D):
+class SymmetricExponential1DSubMesh(SubMesh1D):
     """
     A class to generate a submesh on a 1D domain in which the points are clustered
     close to the boundaries using an exponential formula on the interval [a,b].
     The first half of the interval is meshed using the gridpoints
 
    .. math::
-    x_{k} = (\\frac{b}{2}-a) + \\frac{\\exp{\\alpha k / N} - 1}{\\exp{\\alpha} - 1}} + a,
+    x_{k} = (b/2-a) + \\frac{\\exp{\\alpha k / N} - 1}{\\exp{\\alpha} - 1}} + a,
 
     for k = 1, ..., N, where N is the number of nodes. The grid spacing is then
-    reflected to contruct the grid on the full interval [a,b].
-
+    reflected to contruct the grid on the full interval [a,b]. Here alpha is
+    a stretching factor. As the number of gridpoints tends to infinity, the ratio
+    of the largest and smallest grid cells tends to exp(alpha).
 
     Parameters
     ----------
@@ -135,8 +135,8 @@ def __init__(self, lims, npts):
         spatial_lims = lims[spatial_var]
         npts = npts[spatial_var.id]
 
-        # Strech factor. TODO: allows parameters to be passed to mesh
-        alpha = 1
+        # Strech factor. TODO: allow parameters to be passed to mesh
+        alpha = 2.3 / 2
 
         # Mesh half-interval [a, b/2]
         if npts % 2 == 0:
@@ -162,3 +162,95 @@ def __init__(self, lims, npts):
         coord_sys = spatial_var.coord_sys
 
         super().__init__(edges, coord_sys=coord_sys)
+
+
+class LeftExponential1DSubMesh(SubMesh1D):
+    """
+    A class to generate a submesh on a 1D domain in which the points are clustered
+    close to the left boundary using an exponential formula on the interval [a,b].
+    The gridpoints are given by
+
+   .. math::
+    x_{k} = (b-a) + \\frac{\\exp{\\alpha k / N} - 1}{\\exp{\\alpha} - 1}} + a,
+
+    for k = 1, ..., N, where N is the number of nodes. Here alpha is
+    a stretching factor. As the number of gridpoints tends to infinity, the ratio
+    of the largest and smallest grid cells tends to exp(alpha).
+
+    Parameters
+    ----------
+    lims : dict
+        A dictionary that contains the limits of the spatial variables
+    npts : dict
+        A dictionary that contains the number of points to be used on each
+        spatial variable. Note: the number of nodes (located at the cell centres)
+        is npts, and the number of edges is npts+1.
+    """
+
+    def __init__(self, lims, npts):
+
+        # check that only one variable passed in
+        if len(lims) != 1:
+            raise pybamm.GeometryError("lims should only contain a single variable")
+
+        spatial_var = list(lims.keys())[0]
+        spatial_lims = lims[spatial_var]
+        npts = npts[spatial_var.id]
+
+        # Strech factor. TODO: allow parameters to be passed to mesh
+        alpha = 2.3
+
+        ii = np.array(range(0, npts + 1))
+        a = spatial_lims["min"]
+        b = spatial_lims["max"]
+        edges = (b - a) * (np.exp(alpha * ii / npts) - 1) / (np.exp(alpha) - 1) + a
+
+        coord_sys = spatial_var.coord_sys
+
+        super().__init__(edges, coord_sys=coord_sys)
+
+
+class RightExponential1DSubMesh(SubMesh1D):
+    """
+    A class to generate a submesh on a 1D domain in which the points are clustered
+    close to the right boundary using an exponential formula on the interval [a,b].
+    The gridpoints are given by
+
+   .. math::
+    x_{k} = (b-a) + \\frac{\\exp{-\\alpha k / N} - 1}{\\exp{-\\alpha} - 1}} + a,
+
+    for k = 1, ..., N, where N is the number of nodes. Here alpha is
+    a stretching factor. As the number of gridpoints tends to infinity, the ratio
+    of the largest and smallest grid cells tends to exp(alpha).
+
+    Parameters
+    ----------
+    lims : dict
+        A dictionary that contains the limits of the spatial variables
+    npts : dict
+        A dictionary that contains the number of points to be used on each
+        spatial variable. Note: the number of nodes (located at the cell centres)
+        is npts, and the number of edges is npts+1.
+    """
+
+    def __init__(self, lims, npts):
+
+        # check that only one variable passed in
+        if len(lims) != 1:
+            raise pybamm.GeometryError("lims should only contain a single variable")
+
+        spatial_var = list(lims.keys())[0]
+        spatial_lims = lims[spatial_var]
+        npts = npts[spatial_var.id]
+
+        # Strech factor. TODO: allow parameters to be passed to mesh
+        alpha = 2.3
+
+        ii = np.array(range(0, npts + 1))
+        a = spatial_lims["min"]
+        b = spatial_lims["max"]
+        edges = (b - a) * (np.exp(-alpha * ii / npts) - 1) / (np.exp(-alpha) - 1) + a
+
+        coord_sys = spatial_var.coord_sys
+
+        super().__init__(edges, coord_sys=coord_sys)