set initial particle concentration distribution in DFN

[xyyrly]

Dear team,

I am trying to assign a positive particle concentration distribution as initial c_s_p field. The idea is to interpolant a 2D(x,r) initial c_s_p matrix(which is solution["Positive particle concentration [mol.m-3]"].entries[:,:,-1]) to spatial variables(x_p, r_var). Here x_p is spacial variable in domain"positive electrode" and r_var is spacial variable in domain "positive particle":

xxx = (x[40:,0], r[:,0])
x_p= pybamm.SpatialVariable("x", domain="positive electrode")
r_var = pybamm.SpatialVariable("r", domain="positive particle")
model.initial_conditions[c_s_p]=pybamm.Interpolant(xxx, solution["Positive particle concentration [mol.m-3]"].entries[:,:,0], (x_p, r_var), interpolator="linear")

One error was reported and said x_p and r_var should have the same domain. How should I set x_p? Or are there methods in pybamm to assign an initial distribution？

Thank you so much for your help!

[brosaplanella]

I think there are two main issues here. First, that the order of the variables in the function should be r first and x later. Second, if you are using standard PyBaMM models you can use the standard spatial variables, e.g. define

x_p = pybamm.standard_spatial_vars.x_p

[xyyrly]

Yesss!
Maybe I can take the positive particle concentration part from DFN for test:

import pybamm
import numpy as np

model = pybamm.BaseModel()
param = pybamm.LithiumIonParameters()
T_s_p=273.15+25
c_s_p = pybamm.Variable(
    "Positive particle concentration [mol.m-3]",
    domain="positive particle",
    auxiliary_domains={"secondary": "positive electrode"},
)
c_s_surf_p = pybamm.surf(c_s_p)
N_s_p = -param.p.prim.D(c_s_p, T_s_p) * pybamm.grad(c_s_p)
model.rhs[c_s_p] = -pybamm.div(N_s_p)
model.boundary_conditions[c_s_p] = {
    "left": (pybamm.Scalar(0), "Neumann"),
    "right": (
        -1 / (param.F * param.p.prim.D(c_s_surf_p, T_s_p)),#just for test
        "Neumann",
    ),
}
rr=np.array([1.6250e-07, 4.8750e-07, 8.1250e-07, 1.1375e-06, 1.4625e-06,
       1.7875e-06, 2.1125e-06, 2.4375e-06, 2.7625e-06, 3.0875e-06,
       3.4125e-06, 3.7375e-06, 4.0625e-06, 4.3875e-06, 4.7125e-06,
       5.0375e-06, 5.3625e-06, 5.6875e-06, 6.0125e-06, 6.3375e-06])
xx=np.array([9.5350e-05, 9.8050e-05, 1.0075e-04, 1.0345e-04, 1.0615e-04,
       1.0885e-04, 1.1155e-04, 1.1425e-04, 1.1695e-04, 1.1965e-04,
       1.2235e-04, 1.2505e-04, 1.2775e-04, 1.3045e-04, 1.3315e-04,
       1.3585e-04, 1.3855e-04, 1.4125e-04, 1.4395e-04, 1.4665e-04])
xxx = (rr,xx)
r_p=pybamm.standard_spatial_vars.r_p
x_p = pybamm.standard_spatial_vars.x_p
Initial_c_s_p=12630.8*np.ones((20,20))#just for test, in real case it is not a constant in the whole domain
model.initial_conditions[c_s_p]=pybamm.Interpolant(xxx, Initial_c_s_p, (r_p,x_p), interpolator="linear")

[xyyrly]

Hi brosaplanella,

Thank you for your time. Could this example work?

[brosaplanella]

Yes, that minimum working example is helpful. The issue is with the interpolant, which requires the variables to have the same domains. Not sure if this is an intended feature or a bug. @rtimms @tinosulzer, thoughts?

[valentinsulzer]

Try x_p = pybamm.PrimaryBroadcast(pybamm.standard_spatial_vars.x_p, "positive particle")

[xyyrly]

Thank you again! I tried x_p = pybamm.PrimaryBroadcast(pybamm.standard_spatial_vars.x_p, "positive particle"). Now the domain is correct, but the dimension of x_p and r_p are both (400,1), instead of (20,1). the values of x_p look like:
array([[9.5350e-05],
[9.5350e-05],
[9.5350e-05],
[9.5350e-05],
[9.5350e-05],
[9.5350e-05],
[9.5350e-05],
[9.5350e-05],
[9.5350e-05],
[9.5350e-05],
[9.5350e-05],
[9.5350e-05],
[9.5350e-05],
[9.5350e-05],
[9.5350e-05],
[9.5350e-05],
[9.5350e-05],
[9.5350e-05],
[9.5350e-05],
[9.5350e-05],
[9.8050e-05],
[9.8050e-05],
[9.8050e-05],
[9.8050e-05],
[9.8050e-05],
[9.8050e-05],
[9.8050e-05],
[9.8050e-05],
[9.8050e-05],
[9.8050e-05],
[9.8050e-05],
[9.8050e-05],
[9.8050e-05],
[9.8050e-05],
[9.8050e-05],
[9.8050e-05],
[9.8050e-05],
[9.8050e-05],
[9.8050e-05],
[9.8050e-05],
[1.0075e-04],
[1.0075e-04],
[1.0075e-04],
[1.0075e-04],
[1.0075e-04],
[1.0075e-04],
[1.0075e-04],
[1.0075e-04],
[1.0075e-04],
[1.0075e-04],
[1.0075e-04],
[1.0075e-04],
[1.0075e-04],
[1.0075e-04],
[1.0075e-04],
[1.0075e-04],
[1.0075e-04],
[1.0075e-04],
[1.0075e-04],
[1.0075e-04],
[1.0345e-04],
[1.0345e-04],
[1.0345e-04],
[1.0345e-04],
[1.0345e-04],
[1.0345e-04],
[1.0345e-04],
[1.0345e-04],
[1.0345e-04],
[1.0345e-04],
[1.0345e-04],
[1.0345e-04],
[1.0345e-04],
[1.0345e-04],
[1.0345e-04],
[1.0345e-04],
[1.0345e-04],
[1.0345e-04],
[1.0345e-04],
[1.0345e-04],
[1.0615e-04],
[1.0615e-04],
[1.0615e-04],
[1.0615e-04],
[1.0615e-04],
[1.0615e-04],
[1.0615e-04],
[1.0615e-04],
[1.0615e-04],
[1.0615e-04],
[1.0615e-04],
[1.0615e-04],
[1.0615e-04],
[1.0615e-04],
[1.0615e-04],
[1.0615e-04],
[1.0615e-04],
[1.0615e-04],
[1.0615e-04],
[1.0615e-04],
[1.0885e-04],
[1.0885e-04],
[1.0885e-04],
[1.0885e-04],
[1.0885e-04],
[1.0885e-04],
[1.0885e-04],
[1.0885e-04],
[1.0885e-04],
[1.0885e-04],
[1.0885e-04],
[1.0885e-04],
[1.0885e-04],
[1.0885e-04],
[1.0885e-04],
[1.0885e-04],
[1.0885e-04],
[1.0885e-04],
[1.0885e-04],
[1.0885e-04],
[1.1155e-04],
[1.1155e-04],
[1.1155e-04],
[1.1155e-04],
[1.1155e-04],
[1.1155e-04],
[1.1155e-04],
[1.1155e-04],
[1.1155e-04],
[1.1155e-04],
[1.1155e-04],
[1.1155e-04],
[1.1155e-04],
[1.1155e-04],
[1.1155e-04],
[1.1155e-04],
[1.1155e-04],
[1.1155e-04],
[1.1155e-04],
[1.1155e-04],
[1.1425e-04],
[1.1425e-04],
[1.1425e-04],
[1.1425e-04],
[1.1425e-04],
[1.1425e-04],
[1.1425e-04],
[1.1425e-04],
[1.1425e-04],
[1.1425e-04],
[1.1425e-04],
[1.1425e-04],
[1.1425e-04],
[1.1425e-04],
[1.1425e-04],
[1.1425e-04],
[1.1425e-04],
[1.1425e-04],
[1.1425e-04],
[1.1425e-04],
[1.1695e-04],
[1.1695e-04],
[1.1695e-04],
[1.1695e-04],
[1.1695e-04],
[1.1695e-04],
[1.1695e-04],
[1.1695e-04],
[1.1695e-04],
[1.1695e-04],
[1.1695e-04],
[1.1695e-04],
[1.1695e-04],
[1.1695e-04],
[1.1695e-04],
[1.1695e-04],
[1.1695e-04],
[1.1695e-04],
[1.1695e-04],
[1.1695e-04],
[1.1965e-04],
[1.1965e-04],
[1.1965e-04],
[1.1965e-04],
[1.1965e-04],
[1.1965e-04],
[1.1965e-04],
[1.1965e-04],
[1.1965e-04],
[1.1965e-04],
[1.1965e-04],
[1.1965e-04],
[1.1965e-04],
[1.1965e-04],
[1.1965e-04],
[1.1965e-04],
[1.1965e-04],
[1.1965e-04],
[1.1965e-04],
[1.1965e-04],
[1.2235e-04],
[1.2235e-04],
[1.2235e-04],
[1.2235e-04],
[1.2235e-04],
[1.2235e-04],
[1.2235e-04],
[1.2235e-04],
[1.2235e-04],
[1.2235e-04],
[1.2235e-04],
[1.2235e-04],
[1.2235e-04],
[1.2235e-04],
[1.2235e-04],
[1.2235e-04],
[1.2235e-04],
[1.2235e-04],
[1.2235e-04],
[1.2235e-04],
[1.2505e-04],
[1.2505e-04],
[1.2505e-04],
[1.2505e-04],
[1.2505e-04],
[1.2505e-04],
[1.2505e-04],
[1.2505e-04],
[1.2505e-04],
[1.2505e-04],
[1.2505e-04],
[1.2505e-04],
[1.2505e-04],
[1.2505e-04],
[1.2505e-04],
[1.2505e-04],
[1.2505e-04],
[1.2505e-04],
[1.2505e-04],
[1.2505e-04],
[1.2775e-04],
[1.2775e-04],
[1.2775e-04],
[1.2775e-04],
[1.2775e-04],
[1.2775e-04],
[1.2775e-04],
[1.2775e-04],
[1.2775e-04],
[1.2775e-04],
[1.2775e-04],
[1.2775e-04],
[1.2775e-04],
[1.2775e-04],
[1.2775e-04],
[1.2775e-04],
[1.2775e-04],
[1.2775e-04],
[1.2775e-04],
[1.2775e-04],
[1.3045e-04],
[1.3045e-04],
[1.3045e-04],
[1.3045e-04],
[1.3045e-04],
[1.3045e-04],
[1.3045e-04],
[1.3045e-04],
[1.3045e-04],
[1.3045e-04],
[1.3045e-04],
[1.3045e-04],
[1.3045e-04],
[1.3045e-04],
[1.3045e-04],
[1.3045e-04],
[1.3045e-04],
[1.3045e-04],
[1.3045e-04],
[1.3045e-04],
[1.3315e-04],
[1.3315e-04],
[1.3315e-04],
[1.3315e-04],
[1.3315e-04],
[1.3315e-04],
[1.3315e-04],
[1.3315e-04],
[1.3315e-04],
[1.3315e-04],
[1.3315e-04],
[1.3315e-04],
[1.3315e-04],
[1.3315e-04],
[1.3315e-04],
[1.3315e-04],
[1.3315e-04],
[1.3315e-04],
[1.3315e-04],
[1.3315e-04],
[1.3585e-04],
[1.3585e-04],
[1.3585e-04],
[1.3585e-04],
[1.3585e-04],
[1.3585e-04],
[1.3585e-04],
[1.3585e-04],
[1.3585e-04],
[1.3585e-04],
[1.3585e-04],
[1.3585e-04],
[1.3585e-04],
[1.3585e-04],
[1.3585e-04],
[1.3585e-04],
[1.3585e-04],
[1.3585e-04],
[1.3585e-04],
[1.3585e-04],
[1.3855e-04],
[1.3855e-04],
[1.3855e-04],
[1.3855e-04],
[1.3855e-04],
[1.3855e-04],
[1.3855e-04],
[1.3855e-04],
[1.3855e-04],
[1.3855e-04],
[1.3855e-04],
[1.3855e-04],
[1.3855e-04],
[1.3855e-04],
[1.3855e-04],
[1.3855e-04],
[1.3855e-04],
[1.3855e-04],
[1.3855e-04],
[1.3855e-04],
[1.4125e-04],
[1.4125e-04],
[1.4125e-04],
[1.4125e-04],
[1.4125e-04],
[1.4125e-04],
[1.4125e-04],
[1.4125e-04],
[1.4125e-04],
[1.4125e-04],
[1.4125e-04],
[1.4125e-04],
[1.4125e-04],
[1.4125e-04],
[1.4125e-04],
[1.4125e-04],
[1.4125e-04],
[1.4125e-04],
[1.4125e-04],
[1.4125e-04],
[1.4395e-04],
[1.4395e-04],
[1.4395e-04],
[1.4395e-04],
[1.4395e-04],
[1.4395e-04],
[1.4395e-04],
[1.4395e-04],
[1.4395e-04],
[1.4395e-04],
[1.4395e-04],
[1.4395e-04],
[1.4395e-04],
[1.4395e-04],
[1.4395e-04],
[1.4395e-04],
[1.4395e-04],
[1.4395e-04],
[1.4395e-04],
[1.4395e-04],
[1.4665e-04],
[1.4665e-04],
[1.4665e-04],
[1.4665e-04],
[1.4665e-04],
[1.4665e-04],
[1.4665e-04],
[1.4665e-04],
[1.4665e-04],
[1.4665e-04],
[1.4665e-04],
[1.4665e-04],
[1.4665e-04],
[1.4665e-04],
[1.4665e-04],
[1.4665e-04],
[1.4665e-04],
[1.4665e-04],
[1.4665e-04],
[1.4665e-04]]
So when I do interpolant, it returns an error: The points in dimension 0 must be strictly ascending or descending
Any suggestions for this?