PyBaMM [Bug]:

PyBaMM Version

24.1

Python Version

3.8.19

Describe the bug

I encountered a shape mismatch issue when using the grad function in PyBaMM to calculate the gradient of a variable in a simple model. The specific problem arises due to a mismatch between the shape of the initial conditions and the shape of the rhs after discretization.

Steps to Reproduce

The code I am using is as follows:

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

# Create the model
model = pybamm.BaseModel()

# Define variable u in the x-dimension
u = pybamm.Variable("u", domain=["negative electrode"])

# Define spatial variable x and time variable t
x = pybamm.SpatialVariable("x", domain=["negative electrode"], coord_sys="cartesian")
t = pybamm.t

# Compute the gradient using PyBaMM's grad function
grad_u = pybamm.grad(u)

# Add time evolution equation to the model
model.rhs = {u: grad_u}

# Apply simple Dirichlet boundary conditions
model.boundary_conditions = {
    u: {"left": (pybamm.Scalar(1), "Dirichlet"), "right": (pybamm.Scalar(0), "Dirichlet")}
}

# Set initial conditions, ensuring the shape matches the grid
initial_value_array = np.full(52, 0.5)  # Create an array of 52 elements to match the shape of rhs
initial_conditions = pybamm.Vector(initial_value_array, domain=["negative electrode"])
model.initial_conditions = {u: initial_conditions}

# Register u as an output variable
model.variables = {"u": u}

# Define geometry and mesh with 52 grid points
geometry = {
    "negative electrode": {x: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}}
}
submesh_types = {"negative electrode": pybamm.Uniform1DSubMesh}
var_pts = {x: 52}  # Set 52 grid points to match initial_conditions and rhs
mesh = pybamm.Mesh(geometry, submesh_types, var_pts)

# Create discretization object
spatial_methods = {"negative electrode": pybamm.FiniteVolume()}
disc = pybamm.Discretisation(mesh, spatial_methods)

# Discretize the model
disc.process_model(model)

# Define the time solver and set the time steps
solver = pybamm.ScipySolver()
t_eval = np.linspace(0, 1, 100)  # 100 time steps from 0 to 1 second
solution = solver.solve(model, t_eval)

# Extract and plot the solution
u_sol = solution["u"].entries

# Ensure x and y dimensions match
x_points = mesh["negative electrode"].nodes

# Print the shape of the solution
print("u_values shape:", u_sol.shape)

Relevant log output

pybamm.expression_tree.exceptions.ModelError: rhs and initial conditions must have the same shape after discretisation but rhs.shape = (52, 1) and initial_conditions.shape = (51, 1) for variable 'u'.

Aug 23 '24 07:08 gaoyankk

Typically we solve equations of the type div(grad(u))=0, not grad(u)=0. What is your full system of equations?

Aug 23 '24 23:08 valentinsulzer

Typically we solve equations of the type div(grad(u))=0, not grad(u)=0. What is your full system of equations?

I want to slove this equations . Could you have some good plans ? $$ \frac{\partial \eta}{\partial x} = \frac{-i_s^-}{\sigma^{\text{eff},-}} + \frac{i_{\text{app}} - i_s^-}{\kappa^{\text{eff},-}(c_e)} - \frac{RT}{F}(1 - 2t_+)\frac{\partial \ln(c_e)}{\partial x} - \frac{\partial \text{OCV-}(c_s^-)}{\partial x} $$ pybamm

Aug 27 '24 07:08 gaoyankk

Can you share the full set of equations? There should be more. For example, you have two boundary conditions for u but only a first-order derivative. Is there an equation for div(grad(u))?

Aug 28 '24 04:08 valentinsulzer

Can you share the full set of equations? There should be more. For example, you have two boundary conditions for u but only a first-order derivative. Is there an equation for div(grad(u))?

sorry,I don't have equation for div(grad(u)) . I just want to solv the equation for grad. Do you have a plan to convert grad to div ? A simple example will also do .

Aug 29 '24 06:08 gaoyankk

It looks like there hasn't been a reply in 30 days, so I'm closing this issue.

Nov 22 '24 00:11 github-actions[bot]

I am running into the same problem. Here is a document of the model I want to implement: Here is my code:

import pybamm
import matplotlib
import importlib.util
import numpy as np

TODO: Implement the single well toy problem for a constant (pre-defined) potential 1st (potentiostatic)
TODO: Once above tested and works implement with constant (pre-defined) current (galvanostatic)

# define base model
model = pybamm.BaseModel()

# define parameters
F = pybamm.constants.F
R = pybamm.constants.R
T = pybamm.Parameter("Temperature [K]")
A = pybamm.Parameter("Particle surface area [m2]")
Vol = pybamm.Parameter("Particle volume [m3]")
i_0 = pybamm.Parameter("Reference exchange current density [A.m-2]")
phi = pybamm.Parameter("Potential [V]")  # define as a constant for now
c_max = pybamm.Parameter("Maximum particle stoichiometry")

# define spatial variable
c = pybamm.SpatialVariable("c", domain="positive electrode", coord_sys="cartesian")

# define variable
g = pybamm.Variable("Probability density", domain="positive electrode")

# define dependant equations
# chemical potential
mew = R * T * pybamm.log(c / (c_max - c))
# wind function
v = -((A * i_0) / (F * Vol)) * ((c / c_max) ** 0.5) * ((1 - (c / c_max)) ** 0.5) * pybamm.sinh(
    (1 / (2 * R * T)) * (F * phi + mew)
)
# flux
flux = g * v
dgdt = -pybamm.grad(flux)
model.rhs[g] = dgdt

# boundary conditions
# use dirichlet for both sides for g
model.boundary_conditions[g] = {
    "left": (pybamm.Scalar(0), "Dirichlet"),
    "right": (pybamm.Scalar(0), "Dirichlet"),
}

# initial conditions
# ===========================
# Simple initial condition
# ===========================
model.initial_conditions[g] = pybamm.Scalar(1)

# model variables
model.variables = {
    "Probability density": g,
}

# ===========================
# Model overall set-up
# ===========================
# define parameter values
param = pybamm.ParameterValues(
    {
        "Temperature [K]": 298,
        "Particle surface area [m2]": 5.03e-13,
        "Particle volume [m3]": 3.35e-20,
        "Reference exchange current density [A.m-2]": 8,
        "Potential [V]": 0.35,  # (for potentiostatic condition)
        "Maximum particle stoichiometry": 1,

    }
)

# define geometry
geometry = {
    "positive electrode": {c: {"min": 0, "max": c_max}}
}

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

# mesh
submesh_types = {"positive electrode": pybamm.Uniform1DSubMesh}
var_pts = {c: 20}
mesh = pybamm.Mesh(geometry, submesh_types, var_pts)

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

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

# solve model
solver = pybamm.IDAKLUSolver()
t = np.linspace(0, 1, 100)
solution = solver.solve(model, t)
solution.plot()

I have tried multiple workarounds using pybamm.div with PrimaryBroadcastToEdges as well as the node_to_edge method for the finite volume, but none of those worked. Hence, I am re-opening this issue to ask for any advice. Thanks

Jan 10 '25 14:01 Dharshannan

UPDATE. I have fixed it and found a workaround, here is my code

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

# TODO: Implement the single well toy problem for a constant (pre-defined) potential 1st (potentiostatic)
# TODO: Once above tested and works implement with constant (pre-defined) current (galvanostatic)

# define base model
model = pybamm.BaseModel()

# define parameters
F = pybamm.constants.F
R = pybamm.constants.R
T = pybamm.Parameter("Temperature [K]")
A = pybamm.Parameter("Particle surface area [m2]")
Vol = pybamm.Parameter("Particle volume [m3]")
i_0 = pybamm.Parameter("Reference exchange current density [A.m-2]")
phi = pybamm.Parameter("Potential [V]")  # define as a constant for now
c_max = pybamm.Parameter("Maximum particle stoichiometry")

# define spatial variable
c = pybamm.SpatialVariable(
    "c", domains={"primary": "positive electrode"}, coord_sys="cartesian"
)

# define variable
# g = pybamm.Variable("Probability density", domain="positive electrode")
g = pybamm.Variable(
    "Probability density", domains={"primary": "positive electrode"}
)

# define dependant equations
# chemical potential
mew = R * T * pybamm.log(c / (c_max - c))
# wind function
v = -((A * i_0) / (F * Vol)) * ((c / c_max) ** 0.5) * ((1 - (c / c_max)) ** 0.5) * pybamm.sinh(
    (1 / (2 * R * T)) * (F * phi + mew)
)

# Define the governing equation
# velocity_on_edges = pybamm.PrimaryBroadcastToEdges(v, ["particle"])
flux_on_edges = pybamm.upwind(g) * v * (v > 0) + pybamm.downwind(g) * v * (v < 0)
model.rhs = {g: -pybamm.div(flux_on_edges)}

# boundary conditions
# use dirichlet for both sides for g
model.boundary_conditions[g] = {
    "left": (pybamm.Scalar(0), "Dirichlet"),
    "right": (pybamm.Scalar(0), "Dirichlet"),
}

# initial conditions
# use a gaussian distribution for g initial or (* simpler to just set to 1)
# ===========================
# Gaussian initial condition
# ===========================
A_norm = 1.8305
alpha = 10
c_centre = 0.5
gauss_dist = A_norm * pybamm.exp(-alpha * ((c - c_centre) ** 2))
# model.initial_conditions[g] = gauss_dist

# # ===========================
# # Simple initial condition
# # ===========================
model.initial_conditions[g] = pybamm.Scalar(1)

# model variables
model.variables = {
    "Probability density": g,
}

# ===========================
# Model overall set-up
# ===========================
# define parameter values
param = pybamm.ParameterValues(
    {
        "Temperature [K]": 298,
        "Particle surface area [m2]": 1,
        "Particle volume [m3]": 1,
        "Reference exchange current density [A.m-2]": 1,
        "Potential [V]": -0.025,  # (for potentiostatic condition)
        "Maximum particle stoichiometry": 1,

    }
)

# define geometry
geometry = {
    "positive electrode": {c: {"min": 0, "max": c_max}}
}

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

# mesh
submesh_types = {"positive electrode": pybamm.Uniform1DSubMesh}
var_pts = {c: 200}
mesh = pybamm.Mesh(geometry, submesh_types, var_pts)

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

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

# solve model
solver = pybamm.IDAKLUSolver()
t = np.linspace(0, 1e6, 100)
solution = solver.solve(model, t)
prob_g = np.asarray(solution["Probability density"].entries)

# Select every 30th time interval (you can adjust this if needed)
time_intervals = np.arange(0, prob_g.shape[1], 30)  # Every 30th time step

# Create a plot
plt.figure(figsize=(8, 6))

# Loop over the selected time intervals
for t_idx in time_intervals:
    # Plot the spatial profile of 'g' at the selected time step
    plt.plot(np.linspace(0, 1, prob_g.shape[0]), prob_g[:, t_idx], label=f"t = {t_idx}")

# Add labels and title
plt.xlabel('Space (Normalized)')
plt.ylabel('Probability Density (g)')
plt.title('Time Evolution of Probability Density')
# Add a legend
plt.legend(loc='center left', bbox_to_anchor=(1.01, 0.5), borderaxespad=0)
# Show the plot
plt.show()

Jan 10 '25 18:01 Dharshannan

Maybe should consider writing a jupyter notebook on how to address similar equation implementations within PyBaMM

Jan 10 '25 18:01 Dharshannan