Custom model ShapeError debugging

[isaacbasil]

Hello everyone,

I am trying to program a custom model in pybamm. It simulates a full Li-ion cell, and solves the mass and charge balances on the macroscale (I do not discretise within electrode particles). I have been replicating the basic_dfn file in the pybamm package as much as possible and using the examples on the webpage. However, when it comes to the discretisation step:

disc = pybamm.Discretisation(mesh, spatial_methods)
        `disc.process_model(self.model)

I am getting a ShapeError. From debugging, it looks like a vector of length 61 called disc_left and a MatrixMultiplication of shape (59,1) called disc_right are created by the discretisation script, which I guess is the source of the error, but I don't know how to fix it. Any pointers on things I could try would be great!

import pybamm
import numpy as np
import matplotlib
import matplotlib.pyplot as plt

matplotlib.use('TkAgg')
from params_2eq_model import parameter_values


class create_2eq_model:
    def __init__(self):
        # import parameters from the newman model
        param = parameter_values

        # adding physical params to dictionary
        self.R = pybamm.Parameter('Ideal gas constant [J.K-1.mol-1]')
        self.F = pybamm.Parameter('Faraday constant [C.mol-1]')

        ### Define parameters  ###
        # simulation params
        I = pybamm.Parameter('Current function [A]')
        min_voltage = pybamm.Parameter('Minimum voltage [V]')  # the min voltage at which to cut the simulation [V]
        max_voltage = pybamm.Parameter('Maximum voltage [V]')  # the max voltage at which to cut the simulation [V]

        min_voltage = 2.5  # the min voltage at which to cut the simulation [V]
        max_voltage = 5  # the max voltage at which to cut the simulation [V]

        # global params
        self.T = pybamm.Parameter('Temperature [K]')

        # define negative electrode params
        eps_s_n = pybamm.Parameter('Negative electrode active material volume fraction')
        brugg_s_n = pybamm.Parameter('Negative electrode Bruggeman coefficient (electrode)')
        # tau_s_n = eps_s_n ** (1 - brugg_s_n) #calculate tau from brugg coeff
        tau_s_n = pybamm.Parameter('Negative electrode tortuosity (electrode)')
        D_s_n = pybamm.Parameter('Negative electrode diffusivity [m2.s-1]')
        c_n_max = pybamm.Parameter('Maximum concentration in negative electrode [mol.m-3]')
        sigma_n = pybamm.Parameter('Negative electrode conductivity [S.m-1]')

        # define positive electrode params
        eps_s_p = pybamm.Parameter('Positive electrode active material volume fraction')
        brugg_s_p = pybamm.Parameter('Positive electrode Bruggeman coefficient (electrode)')
        # tau_s_p = eps_s_p ** (1 - brugg_s_p) #calculate tau from brugg coeff
        tau_s_p = pybamm.Parameter('Positive electrode tortuosity (electrode)')
        D_s_p = pybamm.Parameter('Positive electrode diffusivity [m2.s-1]')
        c_p_max = pybamm.Parameter('Maximum concentration in positive electrode [mol.m-3]')
        sigma_p = pybamm.Parameter('Positive electrode conductivity [S.m-1]')

        # define electrolyte parameters
        eps_e_n = pybamm.Parameter('Negative electrode porosity')
        eps_e_p = pybamm.Parameter('Positive electrode porosity')
        brugg_e_n = pybamm.Parameter('Negative electrode Bruggeman coefficient (electrolyte)')
        brugg_e_p = pybamm.Parameter('Positive electrode Bruggeman coefficient (electrolyte)')
        # tau_e_n = eps_e_n ** (1 - brugg_e_n) #calculate tau from brugg coeff
        tau_e_n = pybamm.Parameter('Negative electrode tortuosity (electrolyte)')
        # tau_e_p = eps_e_p ** (1 - brugg_e_p) #calculate tau from brugg coeff
        tau_e_p = pybamm.Parameter('Positive electrode tortuosity (electrolyte)')
        D_e = pybamm.Parameter('Electrolyte diffusivity [m2.s-1]')
        transp_no = pybamm.Parameter('Cation transference number')
        sigma_e = pybamm.Parameter('Electrolyte conductivity [S.m-1]')
        tdf = pybamm.Parameter('Thermodynamic factor')
        theta = - 2 * self.R * self.T * (1 - transp_no) / self.F * tdf * sigma_e

        # define separator parameters
        eps_sep = pybamm.Parameter('Separator porosity')
        brugg_sep = pybamm.Parameter('Separator Bruggeman coefficient (electrolyte)')
        # tau_sep = eps_sep ** (1 - brugg_sep)
        tau_sep = pybamm.Parameter("Separator tortuosity (electrolyte)")

        # define specific surface params which weren't given in import params file
        av_n = pybamm.Parameter('Specific surface area of negative electrode [m2.m-3]')
        av_p = pybamm.Parameter('Specific surface area of positive electrode [m2.m-3]')

        # Initial value params
        c_s_n_0 = pybamm.Parameter('Initial concentration in negative electrode [mol.m-3]')
        c_s_p_0 = pybamm.Parameter('Initial concentration in positive electrode [mol.m-3]')
        c_e_0 = pybamm.Parameter('Initial concentration in electrolyte [mol.m-3]')

        # define geometry params
        L_n = pybamm.Parameter('Negative electrode thickness [m]')
        L_s = pybamm.Parameter('Separator thickness [m]')
        L_p = pybamm.Parameter('Positive electrode thickness [m]')

        # create a base model
        self.model = pybamm.BaseModel()

        ### Define variables ###

        # define negative electrode variables
        c_s_n = pybamm.Variable("Negative electrode average concentration [mol.m-3]", domain="negative electrode")
        phi_s_n = pybamm.Variable("Negative electrode potential [V]", domain="negative electrode")
        c_e_n = pybamm.Variable("Negative electrolyte concentration [mol.m-3]", domain="negative electrode")
        phi_e_n = pybamm.Variable("Negative electrolyte potential [V]", domain="negative electrode")

        # define separator varaiables
        c_e_s = pybamm.Variable("Separator electrolyte concentration [mol.m-3]", domain="separator")
        phi_e_s = pybamm.Variable("Separator electrolyte potential [V]", domain="separator")

        # define positive electrode variables
        c_s_p = pybamm.Variable("Positive electrode average concentration [mol.m-3]", domain="positive electrode")
        phi_s_p = pybamm.Variable("Positive electrode potential [V]", domain="positive electrode")
        c_e_p = pybamm.Variable("Positive electrolyte concentration [mol.m-3]", domain="positive electrode")
        phi_e_p = pybamm.Variable("Positive electrolyte potential [V]", domain="positive electrode")

        # concatenate electrolyte variables and parameters, since it can be described with the same PDEs across all domains
        c_e = pybamm.concatenation(c_e_n, c_e_s, c_e_p)
        phi_e = pybamm.concatenation(phi_e_n, phi_e_s, phi_e_p)
        eps_e_n = pybamm.PrimaryBroadcast(eps_e_n,
                                          "negative electrode") 
        eps_sep = pybamm.PrimaryBroadcast(eps_sep,
                                          "separator") 
        eps_e_p = pybamm.PrimaryBroadcast(eps_e_p,
                                          "positive electrode") 
        tau_e_n = pybamm.PrimaryBroadcast(tau_e_n,
                                          "negative electrode") 
        tau_sep = pybamm.PrimaryBroadcast(tau_sep,
                                          "separator")  
        tau_e_p = pybamm.PrimaryBroadcast(tau_e_p,
                                          "positive electrode")  
        eps_e = pybamm.concatenation(eps_e_n, eps_sep, eps_e_p)
        tau_e = pybamm.concatenation(tau_e_n, tau_sep, tau_e_p)

        ### Define functions ###

        j_n_inputs = {"Negative electrolyte concentration [mol.m-3]": c_e_n,
                      "Negative electrode average concentration [mol.m-3]": c_s_n,
                      'Maximum concentration in negative electrode [mol.m-3]': c_n_max,
                      "Negative electrolyte potential [V]": phi_e_n,
                      "Negative electrode potential [V]": phi_s_n,
                      'Temperature [K]': self.T,
                      }

        j_p_inputs = {"Positive electrolyte concentration [mol.m-3]": c_e_p,
                      "Positive electrode average concentration [mol.m-3]": c_s_p,
                      'Maximum concentration in positive electrode [mol.m-3]': c_p_max,
                      "Positive electrolyte potential [V]": phi_e_p,
                      "Positive electrode potential [V]": phi_s_p,
                      'Temperature [K]': self.T,
                      }

        initial_u0_n_inputs = {"Negative electrode initial SOC": (c_s_n_0 / c_n_max)}
        initial_u0_p_inputs = {"Positive electrode initial SOC": (c_s_p_0 / c_p_max)}

        j_ave_n = pybamm.FunctionParameter("Negative electrode average reaction rate [mol.m-2]", j_n_inputs)
        j_ave_p = pybamm.FunctionParameter("Positive electrode average reaction rate [mol.m-2]", j_p_inputs)
        initial_u0_n = pybamm.FunctionParameter("Negative electrode initial OCP [V]", initial_u0_n_inputs)
        initial_u0_p = pybamm.FunctionParameter("Positive electrode initial OCP [V]", initial_u0_p_inputs)

        ### Calculate reaction rates, flux, and current ###

        # calculate reaction rates
        a_j_n = av_n * j_ave_n  # assuming that surface conc represented by average conc
        j_s = pybamm.PrimaryBroadcast(0, "separator")
        a_j_p = av_p * j_ave_p  # assuming that surface conc represented by average conc
        a_j = pybamm.concatenation(a_j_n, j_s, a_j_p)

        # calculate flux
        N_s_n = - (eps_s_n * D_s_n / tau_s_n) * pybamm.grad(c_s_n)  # negative AM flux
        # N_e_n = - (eps_e_n * D_e / tau_e_n) * pybamm.grad(c_e_n)  # negative electrolyte flux
        N_e = - (eps_e * D_e / tau_e) * pybamm.grad(c_e)  # negative electrolyte flux
        N_s_p = - (eps_s_p * D_s_p / tau_s_p) * pybamm.grad(c_s_p)  # negative AM flux

        # calculate current
        i_n = - eps_s_n * (sigma_n / tau_s_n) * pybamm.grad(phi_s_n)  # negative electrode volume averaged current
        i_e = (eps_e / tau_e) * (- sigma_e * pybamm.grad(phi_e) - theta * pybamm.grad(
            pybamm.log(c_e)))  # electrolyte volume averaged current
        i_p = - eps_s_p * (sigma_p / tau_s_p) * pybamm.grad(phi_s_p)  # negative electrode (macroscale) current

        ### Defind governing equations ###

        # define negative electrode mass balance
        dc_s_n_dt = 1 / eps_s_n * (- pybamm.div(N_s_n) - (
                    1 - transp_no) * a_j_n)  # define the rhs equation, assuming averaged conc represents surface conc.
        self.model.rhs = {c_s_n: dc_s_n_dt}  # add the equation to rhs dictionary

        # define electrolyte mass balance
        dc_e_dt = 1 / eps_e * (-pybamm.div(N_e) + (
                    1 - transp_no) * a_j)  # define the rhs equation, assuming averaged conc represents surface conc.
        self.model.rhs[c_e] = dc_e_dt  # add the equation to rhs dictionary

        # define positive electrode mass balance
        dc_s_p_dt = 1 / eps_s_p * (- pybamm.div(N_s_p) - (
                    1 - transp_no) * a_j_p)  # define the rhs equation, assuming averaged conc represents surface conc.
        self.model.rhs[c_s_p] = dc_s_p_dt  # add the equation to rhs dictionary

        # define negative electrode charge balance
        self.model.algebraic[phi_s_n] = - pybamm.div(i_n) - self.F * a_j_n  # TODO: add preconditioning as in basic dfn

        # define electrolyte charge balance
        self.model.algebraic[phi_e] = - pybamm.div(i_e) + self.F * a_j

        # define positive electrode charge balance
        self.model.algebraic[phi_s_p] = - pybamm.div(i_p) - self.F * a_j_p

        ### define initial conditions ###

        self.model.initial_conditions[c_s_n] = c_s_n_0
        self.model.initial_conditions[c_e] = c_e_0
        self.model.initial_conditions[c_s_p] = c_s_p_0
        # Initial conditions must also be provided for algebraic equations, as an
        # initial guess for a root-finding algorithm which calculates consistent initial
        # conditions
        self.model.initial_conditions[phi_s_n] = pybamm.Scalar(0)
        self.model.initial_conditions[phi_s_p] = initial_u0_p - initial_u0_n
        # initial guess for electrolyte phi as minus negative electrode OCP
        self.model.initial_conditions[phi_e] = - initial_u0_n

        ### define boundary conditions ###

        self.model.boundary_conditions[c_s_n] = {
            "left": (pybamm.Scalar(0), "Neumann"),
            "right": (pybamm.Scalar(0), "Neumann"),
        }

        self.model.boundary_conditions[c_s_p] = {
            "left": (pybamm.Scalar(0), "Neumann"),
            "right": (pybamm.Scalar(0), "Neumann"),
        }

        self.model.boundary_conditions[c_e] = {
            "left": (pybamm.Scalar(0), "Neumann"),
            "right": (pybamm.Scalar(0), "Neumann"),
        }

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

        grad_phi_rhs = - I / sigma_p
        self.model.boundary_conditions[phi_s_p] = {
            "left": (pybamm.Scalar(0), "Neumann"),
            "right": (grad_phi_rhs, "Neumann"),
        }

        self.model.boundary_conditions[phi_e] = {
            "left": (pybamm.Scalar(0), "Neumann"),
            "right": (pybamm.Scalar(0), "Neumann"),
        }

        ### Create variables for visualising the solution ###

        voltage = pybamm.boundary_value(phi_s_p, "right")

        self.model.variables = {
            "Negative electrode average concentration [mol.m-3]": c_s_n,
            "Electrolyte concentration [mol.m-3]": c_e,
            "Positive electrode average concentration [mol.m-3]": c_s_p,
            "Current [A]": I,
            "Negative electrode potential [V]": phi_s_n,
            "Electrolyte potential [V]": phi_e,
            "Positive electrode potential [V]": phi_s_p,
            "Voltage [V]": voltage,
            "Time [s]": pybamm.t,
        }

        ### Termination events ###
        # Events specify points at which a solution should terminate
        self.model.events += [
            pybamm.Event("Minimum voltage [V]", voltage - min_voltage),
            pybamm.Event("Maximum voltage [V]", max_voltage - voltage),
        ]

        ### Define geometry and discretise ###

        # define geometry
        x_n = pybamm.SpatialVariable("x_n", domain=["negative electrode"], coord_sys="cartesian")
        x_s = pybamm.SpatialVariable("x_s", domain=["separator"], coord_sys="cartesian")
        x_p = pybamm.SpatialVariable("x_p", domain=["positive electrode"], coord_sys="cartesian")

        geometry = {"negative electrode": {x_n: {"min": pybamm.Scalar(0), "max": L_n}},
                    "separator": {x_s: {"min": L_n, "max": (L_n + L_s)}},
                    "positive electrode": {x_p: {"min": (L_n + L_s), "max": (L_n + L_s + L_p)}},
                    }

        param.process_model(self.model)
        param.process_geometry(geometry)

        # mesh and discretise
        submesh_types = {"negative electrode": pybamm.Uniform1DSubMesh,
                         "separator": pybamm.Uniform1DSubMesh,
                         "positive electrode": pybamm.Uniform1DSubMesh,
                         }

        var_pts = {x_n: 20,
                   x_s: 20,
                   x_p: 20,
                   }

        mesh = pybamm.Mesh(geometry, submesh_types, var_pts)

        # define spatial methods
        spatial_methods = {"negative electrode": pybamm.FiniteVolume(),
                           "separator": pybamm.FiniteVolume(),
                           "positive electrode": pybamm.FiniteVolume(),
                           }

        disc = pybamm.Discretisation(mesh, spatial_methods)
        disc.process_model(self.model)

        # solve
        solver = pybamm.CasadiSolver()
        t = np.linspace(0, 1000, 100)
        solution = solver.solve(self.model, t)


simple_2eq = create_2eq_model()

if you have read this far.. sorry for the code style, I just wanted to get everything working before I started formatting the class :)

[rtimms]

This kind of off-by-two shape error is usually because you are trying to add together a gradient (which is defined on cell edges in the finite volume method) to a variable (which is defined on cell centres), so I'd look at checking your equations to see if you are accidentally doing this anywhere. You can add pybamm.set_logging_level("DEBUG") to the start of your script to get more info (e.g. it will tell you which variable equation it is trying to discretise when you get the error).

[rtimms]

It may also be that you haven't provided boundary conditions for a quantity inside a divergence operator.

[isaacbasil]

Thanks that's really helpful, I am indeed trying to add a gradient and a variable together. To get around that, I suppose I will need to use some kind of broadcast to define the cell center values on the boundaries (or vice versa)?

[isaacbasil]

To follow up on this - the problem seemed to go away when I removed a logarithm from one of the equations. In the expression for i_e, I had a grad(log(c_e)) term. Instead I wrote this as 1/c_e * grad(c_e). After that, things seem to run smoothly. Not sure why that is :)

[valentinsulzer]

PyBaMM looks for boundary conditions with keys matching what's inside grad. In your case you likely had a boundary condition for c_e, not for log(c_e).