Documentation on the integral operator

[isaacbasil]

Hi,

As per my previous issue, I want to start contributing to Pybamm, but do not have any experience contributing to open-source projects so I am looking for an easy entry point.

In my work I have been using pybamm's integral operator. I find the documentation a bit misleading here. Currently it describes the integral operator as follows,

pybamm.Integral(function, integration_variable)

$I = \int_{a}^{b}\!f(u)\,du$

"where $a$ and $b$ are the left-hand and right-hand boundaries of the domain respectively, and $u\in\text{domain}$."

However, I don't think that the differential element actually used in calculating the integral is necessarily du (where u is the integration variable). For example,

$\int_{0}^{R}\!1\,dr = R$

But,

pybamm.Integral(1, r) $\neq R$ (at least in certain cases, as shown below)

It looks like the integral operator actually performs,

$I = \int_{u_{min}}^{u_{max}}\!f(u)\,dq$

Where $u$ is some spatial variable, $u_{min}$ and $u_{max}$ are set by the geometry of the spatial variable, and $dq$ is given by,

$dq = du$ for cartesian coordinates.
$dq = 2\pi u du$ for cylindrical coordinates.
$dq = 4\pi u^2 du$ for spherical coordinates.

I show this in the code below for the cylindrical case, where,

$\int_{0}^{R}\!1\,dq = \pi R^2$

In the code, if you change coord_sys="cylindrical polar" to coord_sys="cartesian", it will give,

$\int_{0}^{R}\!1\,dq = R $

and if you change it to coord_sys="spherical polar" it will be

$\int_{0}^{R}\!1\,dq = \frac{4}{3} \pi r^3 $

As mentioned above, I'm looking for a contribution to make myself - I am just starting an issue first to follow the contribution guidelines. Please let me know your thoughts on this suggestion.

Steps to Reproduce

import pybamm
import numpy as np

model = pybamm.BaseModel()

r_p = pybamm.SpatialVariable("r", domain="positive particle", coord_sys="cylindrical polar")

R_p = 2

one = pybamm.PrimaryBroadcast(1, "positive particle")
one_int = pybamm.Integral(one, r_p)


geometry = {"positive particle": {r_p: {"min": 0, "max": R_p }}}

var_pts = {r_p: 50}

# the PDE describing c is irrelevant, it's just to run the solver 

c = pybamm.Variable("c", domain="positive particle")
D = pybamm.Scalar(1)
dcdt = pybamm.div(D * pybamm.grad(c))
model.rhs[c] = dcdt

model.boundary_conditions = {
c: {"left": (0, "Neumann"),
"right": (1, "Neumann")}, }
model.initial_conditions = {c: 1}

model.variables = {"c": c, "Integral of 1": one_int}

submesh_types = {"positive particle": pybamm.Uniform1DSubMesh}
mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
spatial_methods = {"positive particle": pybamm.FiniteVolume()}
disc = pybamm.Discretisation(mesh, spatial_methods)
disc.process_model(model)
sim = pybamm.Simulation(model)
sim.solve([0, 1])

integral_result = sim.solution["Integral of 1"].entries[0]

pi_r_squared = np.pi * R_p **2
v_sphere =  4 / 3 * np.pi * R_p **3

print("The integral calculated by pybamm is ", integral_result)
print("The value of R is ", R_p)
print("The value of pi * R^2 is ", pi_r_squared)
print("The value of 4/3 * pi * R^3 is ", v_sphere)

[martinjrobins]

thanks @isaacbasil! It would indeed be helpful to provided this additional detail in the docstring wrt the different coordinate systems. Please do make a PR for this