Restructure "Creating Models" notebooks

[brosaplanella]

After reviewing the "Getting Started" notebooks, the next step is to restructure the "Creating Models" ones. The idea is to have a cohesive set of notebooks that build progressively more complex models up to DFN (potentially with degradation). The inspiration is the "12 steps to Navier Stokes" tutorial.

We will not have 12 steps, at least initially, but so far the structure would be:

1. Reservoir model (ODEs)
2. SPM (parabolic PDEs)
3. SPMe (parabolic PDEs in multiple domains, multiscale formulation appears)
4. DFN (parabolic & elliptic PDEs in multiple domains)
5. DFN + degradation (maybe SEI? Includes all of the above, but more complicated)

This is subject to change, as the notebooks are being edited/written (e.g. some might need to be split into two).

This issue would go hand-in-hand with this: https://github.com/pybamm-team/PyBaMM/issues/3839#issuecomment-1966614301

[valentinsulzer]

Do we want all the models to be battery models or should we start with generic models? Pro of battery models is it's more interesting, con is it requires loading in a lot of parameters

If we can have generic models, before 1 could be solving an exponential decay equation (introduce the concept of a variable, rhs, initial conditions, etc) and between 1 and 2 could be diffusion in a cartesian domain.

Other comments:

• Full SPMe is tricky because it requires a complex expression for the voltage
• It might be helpful to include some intermediate models between SPMe and DFN that introduce the coupling of potential as an algebraic equation

[valentinsulzer]

Here are some models I wrote a while back, they don't work any more but might be useful intermediates to explain some of the model structure

import pybamm

class BaseModel(pybamm.lithium_ion.BaseModel):
    def __init__(self, name):
        super().__init__({}, name)
        self.variables = {
            "Time": pybamm.t,
            "x": pybamm.standard_spatial_vars.x,
            "x_n": pybamm.standard_spatial_vars.x_n,
            "x_p": pybamm.standard_spatial_vars.x_p,
            "r_n": pybamm.SpatialVariable("r_n", domain=["negative particle"]),
            "r_p": pybamm.SpatialVariable("r_p", domain=["positive particle"]),
        }

    def new_empty_copy(self):
        return pybamm.BaseModel.new_empty_copy(self)

def spm():
    model = pybamm.lithium_ion.SPM(name="SPM")
    model.variables.update(
        {
            "r_n": pybamm.SpatialVariable("r_n", domain=["negative particle"]),
            "r_p": pybamm.SpatialVariable("r_p", domain=["positive particle"]),
        }
    )
    variables = [
        "X-averaged negative particle concentration",
        "X-averaged positive particle concentration",
    ]
    return model, variables

def spm_no_r():
    model = pybamm.lithium_ion.SPM({"particle": "uniform profile"}, name="SPM_no_r")
    variables = [
        "Average negative particle concentration",
        "Average positive particle concentration",
    ]
    return model, variables

def spme():
    model = pybamm.lithium_ion.SPMe(name="SPMe")
    model.variables.update(
        {
            "r_n": pybamm.SpatialVariable("r_n", domain=["negative particle"]),
            "r_p": pybamm.SpatialVariable("r_p", domain=["positive particle"]),
        }
    )
    variables = [
        "X-averaged negative particle concentration",
        "X-averaged positive particle concentration",
        "Electrolyte concentration",
    ]
    return model, variables

def dfn_no_r():
    model = pybamm.lithium_ion.DFN({"particle": "uniform profile"}, name="DFN_no_r")
    variables = [
        "R-averaged negative particle concentration",
        "R-averaged positive particle concentration",
        "Electrolyte concentration",
        "Electrolyte potential",
        "Negative electrode potential",
        "Positive electrode potential",
    ]
    return model, variables

def dfn():
    model = pybamm.lithium_ion.DFN(name="DFN")
    variables = [
        "Negative particle concentration",
        "Positive particle concentration",
        "Electrolyte concentration",
        "Electrolyte potential",
        "Negative electrode potential",
        "Positive electrode potential",
    ]
    return model, variables

def reduced_c():
    class Model(BaseModel):
        def __init__(self, name="reduced_c"):
            super().__init__(name)
            param = self.param
            c_e = pybamm.Variable(
                "Electrolyte concentration",
                ["negative electrode", "separator", "positive electrode"],
            )

            a_j = pybamm.concatenation(
                pybamm.PrimaryBroadcast(1000 / param.n.L, "negative electrode"),
                pybamm.PrimaryBroadcast(0, "separator"),
                pybamm.PrimaryBroadcast(-1000 / param.p.L, "positive electrode"),
            )
            self.rhs[c_e] = pybamm.div(pybamm.grad(c_e)) + a_j

            self.initial_conditions[c_e] = pybamm.Scalar(1)
            self.boundary_conditions[c_e] = {
                "left": (pybamm.Scalar(0), "Neumann"),
                "right": (pybamm.Scalar(0), "Neumann"),
            }
            self.variables.update({"Electrolyte concentration": c_e})

    model = Model()
    variables = list(model.variables.keys())
    return model, variables

def reduced_c_phi():
    class Model(BaseModel):
        def __init__(self, name="reduced_c_phi"):
            super().__init__(name)
            param = self.param
            c_e = pybamm.Variable(
                "Electrolyte concentration",
                ["negative electrode", "separator", "positive electrode"],
            )
            phi_e = pybamm.Variable(
                "Electrolyte potential",
                ["negative electrode", "separator", "positive electrode"],
            )

            a_j = pybamm.concatenation(
                pybamm.PrimaryBroadcast(1 / param.n.L, "negative electrode"),
                pybamm.PrimaryBroadcast(0, "separator"),
                pybamm.PrimaryBroadcast(-1 / param.p.L, "positive electrode"),
            )
            self.rhs[c_e] = pybamm.div(pybamm.grad(c_e)) + a_j
            self.algebraic[phi_e] = (
                pybamm.div(pybamm.grad(c_e) / c_e - pybamm.grad(phi_e)) - a_j
            )

            self.initial_conditions[c_e] = pybamm.Scalar(1)
            self.initial_conditions[phi_e] = pybamm.Scalar(0)
            self.boundary_conditions[c_e] = {
                "left": (pybamm.Scalar(0), "Neumann"),
                "right": (pybamm.Scalar(0), "Neumann"),
            }
            self.boundary_conditions[phi_e] = {
                "left": (pybamm.Scalar(0), "Dirichlet"),
                "right": (pybamm.Scalar(0), "Neumann"),
            }
            self.variables.update(
                {"Electrolyte concentration": c_e, "Electrolyte potential": phi_e}
            )

    model = Model()
    variables = list(model.variables.keys())
    return model, variables

def reduced_c_phi_j():
    class Model(BaseModel):
        def __init__(self, name="reduced_c_phi_j"):
            super().__init__(name)
            param = self.param
            c_e_n = pybamm.Variable(
                "Negative electrolyte concentration", "negative electrode"
            )
            c_e_s = pybamm.Variable("Separator electrolyte concentration", "separator")
            c_e_p = pybamm.Variable(
                "Positive electrolyte concentration", "positive electrode"
            )
            c_e = pybamm.concatenation(c_e_n, c_e_s, c_e_p)

            phi_e_n = pybamm.Variable(
                "Negative electrolyte potential", "negative electrode"
            )
            phi_e_s = pybamm.Variable("Separator electrolyte potential", "separator")
            phi_e_p = pybamm.Variable(
                "Positive electrolyte potential", "positive electrode"
            )
            phi_e = pybamm.concatenation(phi_e_n, phi_e_s, phi_e_p)

            phi_s_n = pybamm.Variable(
                "Negative electrode potential", "negative electrode"
            )
            phi_s_p = pybamm.Variable(
                "Positive electrode potential", "positive electrode"
            )

            a_j_n = c_e_n**0.5 * pybamm.sinh(phi_s_n - phi_e_n + 1)
            a_j_p = c_e_p**0.5 * pybamm.sinh(phi_s_p - phi_e_p - 4)
            a_j = pybamm.concatenation(
                a_j_n, pybamm.PrimaryBroadcast(0, "separator"), a_j_p
            )

            self.rhs[c_e] = pybamm.div(pybamm.grad(c_e)) + a_j
            self.algebraic[phi_e] = (
                pybamm.div(pybamm.grad(c_e) / c_e - pybamm.grad(phi_e)) - a_j
            )
            self.algebraic[phi_s_n] = pybamm.div(pybamm.grad(phi_s_n)) + a_j_n
            self.algebraic[phi_s_p] = pybamm.div(pybamm.grad(phi_s_p)) + a_j_p

            self.initial_conditions[c_e] = pybamm.Scalar(1)
            self.initial_conditions[phi_e] = pybamm.Scalar(1)
            self.initial_conditions[phi_s_n] = pybamm.Scalar(0)
            self.initial_conditions[phi_s_p] = pybamm.Scalar(4)

            self.boundary_conditions[c_e] = {
                "left": (pybamm.Scalar(0), "Neumann"),
                "right": (pybamm.Scalar(0), "Neumann"),
            }
            self.boundary_conditions[phi_e] = {
                "left": (pybamm.Scalar(0), "Neumann"),
                "right": (pybamm.Scalar(0), "Neumann"),
            }
            self.boundary_conditions[phi_s_n] = {
                "left": (pybamm.Scalar(0), "Dirichlet"),
                "right": (pybamm.Scalar(0), "Neumann"),
            }
            self.boundary_conditions[phi_s_p] = {
                "left": (pybamm.Scalar(0), "Neumann"),
                "right": (pybamm.Scalar(1), "Neumann"),
            }

            self.variables.update(
                {
                    "Electrolyte concentration": c_e,
                    "Electrolyte potential": phi_e,
                    "Negative electrode potential": phi_s_n,
                    "Positive electrode potential": phi_s_p,
                    "Interfacial current density": a_j,
                }
            )

    model = Model()
    variables = list(model.variables.keys())
    return model, variables

def solve_plot(model, variables):
    sim = pybamm.Simulation(model)
    sim.solve([0, 3600])
    sim.plot(variables)

# solve_plot(*spm())
# solve_plot(*spm_no_r())
# solve_plot(*reduced_c())
# solve_plot(*spme())
# solve_plot(*reduced_c_phi())
# solve_plot(*reduced_c_phi_j())
# solve_plot(*dfn_no_r())
# solve_plot(*dfn())