setting ["Ambient temperature [K]"] as pybamm.Interpolant and getting broadcast error #3568

swhan22 · 2023-11-27T18:02:20Z

swhan22
Nov 27, 2023

Hi, I am trying to run a drive cycle simulation. I am able to set ["Current function [A]"] as pybamm.Interpolant and run DFN model without any error. However, as soon as I set ["Ambient temperature [K]"] as Interpolant, I am getting broadcast error (please see below). I am fairly new to python. Could any one give me clues to what could be the issue? Is there any other method to give temperature profile as an input to DFN model? Thank you.🙏

error message:
Cannot find shape (original error: operands could not be broadcast together with shapes (115,1) (41,1) );

I am using the following 1D mesh size:
var_pts = {var.x_n: 50, var.x_s: 15, var.x_p: 50, var.r_n: 20, var.r_p: 20}

Answered by rtimms

Nov 28, 2023

Ambient temperature can now depend on y,z,t, so you'll need to wrap your interpolant function. Somthing like

def T_function(y,z,t):
    return T_interp_function(t)

View full answer

valentinsulzer · 2023-11-27T19:46:33Z

valentinsulzer
Nov 27, 2023
Maintainer

Could you post a minimum reproducible example please?

4 replies

swhan22 Nov 27, 2023
Author

Hi Valentin, Thanks so much for your response. Below is a minimum reproducible example.

As soon as I uncomment line -
para_set["Ambient temperature [K]"] = T_interp_function

I get the broadcast error.

x100 = 0.85
x0 = 0.14595

y100 = 0.18
y0 = 0.88178

c1maxpos = 4900 / 97.2795 * 1e3
c1maxneg = 2250 / 72.066 * 1e3

def Upos(sto):
    
    p1 = -5.9309
    p2 = -1.3492
    p3 = 61.55
    p4 = -5.5573
    p5 = -0.40504
    p6 = -386.94
    p7 = 696.7
    p8 = -509.27
    p9 = 183.82
    p10 = -33.499
    p11 = 7.5365
    p12 = -0.87638
    p13 = 35.559
    p14 = 39.962
    
    u_eq = (
       p1*sto**10 + p2*sto**9 + p3*sto**8 + p4*sto**7 + p5*sto**6 + p6*sto**5
       + p7*sto**4 + p8*sto**3 + p9*sto**2 + p10*sto + p11 + p12*pb.exp(p13*sto**p14)
    )
      
    return u_eq



def Uneg(sto):
    
    p1 = 0.15088
    p2 = 0.5789
    p3 = -60.535
    p4 = 0.48031
    p5 = -521.38
    p6 = 41.674
    p7 = -37.204
    p8 = -0.030289
    p9 = 28.483
    p10 = -3.9778
    p11 = -0.012591
    p12 = 30.114
    p13 = -13.481    
    
    u_eq = (
    p1 + p2*pb.exp(p3*sto) + p4*pb.exp(p5*sto) 
          - pb.exp(p6*sto + p7) + p8*pb.arctan(p9*sto + p10)
          + p11*pb.arctan(p12*sto + p13)
    )

    return u_eq

  
def D1neg(sto, T):
    
    Dneg_constant = 3e-14
    Ea_Dneg = 20e3
    Tref = 40 + 273.15
    
    arrhenius = pb.exp(Ea_Dneg / pb.constants.R * (1 / Tref - 1 / T))
    
    return arrhenius * Dneg_constant


def D1pos(sto, T):
    Dpos_constant = 8e-14
    Ea_Dpos = 20e3
    Tref = 40 + 273.15

    arrhenius = pb.exp(Ea_Dpos / pb.constants.R * (1 / Tref - 1 / T))

    return arrhenius * Dpos_constant


def i0neg(c_e, c_s_surf, c_s_max, T):
    
    kneg = 3e-9
    Ea_kneg = 45e3
    c_e_0 = 1200
    Tref = 40 + 273.15

    arrhenius = pb.exp(Ea_kneg / pb.constants.R * (1/Tref - 1/T))

    return pb.constants.F * kneg * arrhenius * (c_e / c_e_0)**0.5 * c_s_surf**0.5 * (c_s_max - c_s_surf)**0.5

def i0pos(c_e, c_s_surf, c_s_max, T):
    kpos = 4e-9
    Ea_kpos = 45e3
    c_e_0 = 1200
    Tref = 40 + 273.15

    arrhenius = pb.exp(Ea_kpos / pb.constants.R * (1/Tref - 1/T))

    return pb.constants.F * kpos * arrhenius * (c_e / c_e_0)**0.5 * c_s_surf**0.5 * (c_s_max - c_s_surf)**0.5    


def electrolyte_diffusivity(c_e, T):
    
    D00 = 4.43
    D01 = 54
    D10 = 0.22
    Tg0 = 229
    Tg1 = 5
    # mol/m3 to molar
    c_e = c_e / 1000

    # from Valoen w/ D11 = 0
    T_g = Tg0 + Tg1 * c_e
    D_0 = -D00 - D01 / (T - T_g)
    D_1 = -D10

    return (10 ** (D_0 + D_1 * c_e)) * 1e-4

def electrolyte_conductivity(c_e, T):
    # mol/m3 to molar
    c_e = c_e / 1000
    # mS/cm to S/m
    return (1e-3 / 1e-2) * (
        c_e
        * (
            (-10.5 + 0.0740 * T - 6.96e-5 * T**2)
            + c_e * (0.668 - 0.0178 * T + 2.80e-5 * T**2)
            + c_e**2 * (0.494 - 8.86e-4 * T)
        )
        ** 2
    )
    
    
def electrolyte_activity(c_e, T):
    t_plus = 0.36
    Tref = 40 + 273.15

    dlnf_dlnc = (((
        0.601 - 0.24 * (c_e / 1000) ** 0.5
        + 0.982 * (1 - 0.0052 * (T - Tref)) * ((c_e / 1000) ** 1.5)
        ) / (1 - t_plus) - 1)) + 1
    
    return dlnf_dlnc


para = {
    "Contact resistance [Ohm]": 0.1e-3,     
    "Electrode height [m]" : 286e-3,
    "Electrode width [m]" : 72.9e-3,
    "Cell volume [m3]" : 286e-3*72.9e-3*38*2*(60.5 + 91.5 + 15)*1e-6*1,
    # number of electrode pairs/layers     
    "Number of electrodes connected in parallel to make a cell" : 38 * 2,
    "Number of cells connected in series to make a battery" : 1.0,

    # Negative domain
    "Negative electrode OCP [V]" : Uneg,
    "Maximum concentration in negative electrode [mol.m-3]" : c1maxneg,
    "Initial concentration in negative electrode [mol.m-3]" : (x0 + 1 * (x100 - x0)) * c1maxneg,
    "Negative electrode diffusivity [m2.s-1]" : D1neg, # CHECK
    "Negative current collector thickness [m]" : 1.2e-15,
    "Negative electrode thickness [m]" : 91.5e-6,
    "Negative electrode conductivity [S.m-1]" : 100,
    "Negative current collector conductivity [S.m-1]" : 58411000.0*1000,
    "Negative current collector thermal conductivity [W.m-1.K-1]" : 401.0*1000,
    "Negative current collector density [kg.m-3]" : 8960.0,
    "Negative current collector specific heat capacity [J.kg-1.K-1]" : 385.0*0,
    "Negative electrode porosity" : 0.319,
    "Negative electrode active material volume fraction" : 1 - 0.319 -0.03,
    "Negative particle radius [m]" : 6e-06,
    "Negative electrode Bruggeman coefficient (electrolyte)" : 1-pb.log(6)/pb.log(0.319),
    "Negative electrode Bruggeman coefficient (electrode)" : 1, 
    "Negative electrode cation signed stoichiometry" : -1.0,
    "Negative electrode electrons in reaction" : 1.0,
    "Negative electrode charge transfer coefficient" : 0.5,
    "Negative electrode double-layer capacity [F.m-2]" : 0.0,
    "Negative electrode exchange-current density [A.m-2]" : i0neg,
    "Negative electrode density [kg.m-3]" : 2250,
    "Negative electrode specific heat capacity [J.kg-1.K-1]" : 700.0*0 + 1000,
    "Negative electrode thermal conductivity [W.m-1.K-1]" : 1.7*1000,
    "Negative electrode OCP entropic change [V.K-1]" : 0.0,
    
    # Sep domain
    "Separator thickness [m]" : 15e-6,
    "Separator porosity" : 0.45, # CHECK
    "Separator Bruggeman coefficient (electrolyte)" : 1-pb.log(5)/pb.log(0.45),
    "Separator density [kg.m-3]" : 397.0,
    "Separator specific heat capacity [J.kg-1.K-1]" : 700.0*0 + 1000,
    "Separator thermal conductivity [W.m-1.K-1]" : 0.16*1000,
    
    # Pos domain
    "Positive electrode OCP [V]" : Upos,
    "Maximum concentration in positive electrode [mol.m-3]" : c1maxpos, # CHECK
    "Initial concentration in positive electrode [mol.m-3]" : (y0 + 1 * (y100 - y0)) * c1maxpos,
    "Positive electrode diffusivity [m2.s-1]" : D1pos,
    "Positive electrode conductivity [S.m-1]" : 3,
    "Positive electrode porosity" : 0.292, 
    "Positive electrode active material volume fraction" : 1 - 0.292 - 0.03, # CHECK
    "Positive electrode Bruggeman coefficient (electrolyte)" : 1 - pb.log(5)/pb.log(0.292), 
    "Positive particle radius [m]" : 5e-06,
    
    "Positive electrode Bruggeman coefficient (electrode)" : 1, 
    "Positive electrode cation signed stoichiometry" : -1.0,
    "Positive electrode electrons in reaction" : 1.0,
    "Positive electrode charge transfer coefficient" : 0.5,
    "Positive electrode double-layer capacity [F.m-2]" : 0.0,
    "Positive electrode exchange-current density [A.m-2]" : i0pos,
    "Positive electrode density [kg.m-3]" : 4900.0,
    "Positive electrode specific heat capacity [J.kg-1.K-1]" : 700.0*0 + 1080.060,
    "Positive electrode thermal conductivity [W.m-1.K-1]" : 2.1*1000,
    "Positive electrode OCP entropic change [V.K-1]" : 0.0,
    "Positive electrode thickness [m]" : 60.5e-6, # CHECK
    "Positive current collector thickness [m]" : 1.6e-15,
    "Positive current collector conductivity [S.m-1]" : 36914000.0*1000,
    "Positive current collector density [kg.m-3]" : 2700.0,
    "Positive current collector specific heat capacity [J.kg-1.K-1]" : 897.0*0,
    "Positive current collector thermal conductivity [W.m-1.K-1]" : 237.0*1000,
    
    # thermal props
    "Reference temperature [K]" : 40 + 273.15,    
    "Total heat transfer coefficient [W.m-2.K-1]" : 75,
    "Ambient temperature [K]" : 40 + 273.15,
    # "Ambient temperature [K]" : T_interp_function,
    "Cell cooling surface area [m2]" : 286e-3*72.9e-3*2,
    "Cell thermal expansion coefficient [m.K-1]" : 1.1e-06*1,
    "tinterp" : 1,
    "Tinterp" : 1,
        

    # electrolyte
    "Typical electrolyte concentration [mol.m-3]" : 1200.0, 
    "Initial concentration in electrolyte [mol.m-3]" : 1200.0, 
    "Cation transference number" : 0.36, 
    "Thermodynamic factor" : electrolyte_activity,    
    
   
    "Electrolyte conductivity [S.m-1]" : electrolyte_conductivity,
    "Electrolyte diffusivity [m2.s-1]" : electrolyte_diffusivity,    
    
    "para_soc00" : 1,
    "para_D00" : 4.43,
    "para_D01" : 54,    
    "para_D10" : 0.22,
    "para_Tg0" : 228,
    "para_Tg1" : 5,
    
    "para_sf_activity" : 1,
    "para_sf_a1_activity" : 1,

    "para_Ea_D2" : 18e3,
    "para_tauPos" : 5,
    "para_tauNeg" : 6,

    "para_Ea_kpos" : 45e3,
    "para_kpos" : 4e-9,
    "para_sf_kappa" : 1,
    
    "para_Ea_kneg" : 45e3,
    "para_kneg" : 3e-9,
    
    "para_Ea_Dpos" : 20e3,
    "para_Ea_Dneg" : 20e3,

    "para_sf_D1pos" : 1,
    "para_sf_D1neg" : 1,
    "lsthermal" : 0,
    
    "Lower voltage cut-off [V]" : 2.7,
    "Upper voltage cut-off [V]" : 4.4,
    "Initial temperature [K]" : 40 + 273.15,
        
    "Nominal cell capacity [A.h]" : 58,
    
    "Typical current [A]" : 58*3,
    
    "Current function [A]" : 58 * 3,    
    
    "Power function [W]": 100,
    
}

pb.settings.tolerances["D_e__c_e"] = 1e-10
pb.settings.tolerances["kappa_e__c_e"] = 1e-10
pb.settings.tolerances["chi__c_e"] = 1e-10
pb.settings.tolerances["U__c_s"] = 1e-10
pb.settings.tolerances["j0__c_e"] = 1e-10
pb.settings.tolerances["j0__c_s"] = 1e-10
pb.settings.tolerances["macinnes__c_e"] = 1e-10


# drive cycle snippet
dc = np.array([[0, 0,	4.178255081, 39.92520905],
                [0.139264, -184.5974121,	4.036471844, 39.8711853],
                [0.237568,	-184.8833618,	4.03249073,	39.8711853],
                [0.335872,	-185.0022583,	4.029868603,	39.8711853],
                [0.442368,	-185.1472168,	4.027357101,	39.8671608],
                [0.53248,	-185.2461548,	4.025063038,	39.8671608],
                [0.638976,	-185.319519,    4.023223877,	39.8671608],
                [0.73728,	-185.4052734,	4.021364689,	39.8671608],
                [0.843776,	-185.4545288,	4.019886494,	39.8671608],
                [0.94208,	-185.5421753,	4.018260956,	39.8671608],
                [1.040384,	-185.6124878,	4.016605854,	39.8671608],
                [1.138688,	-185.6817627,	4.015146732,	40.02028656],
                [1.236992,	-185.7332153,	4.013901234,	40.02028656],
                [1.343488,	-185.8170166,	4.012485027,	40.02028656],
                [1.441792,	-185.8728027,	4.011345387,	40.02028656],
                [1.540096,	-185.914856,    4.010089397,	40.02028656],
                [1.646592,	-185.9702148,	4.008943081,	40.02028656],
                [1.744896,	-186.0526123,	4.007782459,	40.02028656],
                [1.8432,	-186.1102295,	4.00637722,	    40.02028656],
                [1.949696,	-186.1390991,	4.005376816,	39.96378326],
                [2.048,	    -186.154541,	4.004361153,	39.96378326],
                [2.154496,	-186.2084961,	4.003252983,	39.96378326],
                [2.2528,	-186.2745972,	4.001954556,	39.96378326],
                [2.367488,	-186.3257446,	4.000907421,	39.96378326],
                [2.490368,	-186.373291,	3.99977231,	    39.96378326],
                [2.613248,	-186.4265747,	3.998755217,	39.99684906],
                [2.711552,	-186.4959717,	3.997748375,	39.99684906],
                [2.842624,	-186.5219727,	3.996665239,	39.99684906],
                [2.94912,	-186.5685425,	3.995596886,	39.99684906],
                [3.080192,	-186.6502075,	3.994441509,	39.99684906],
                [3.219456,	-186.6798706,	3.993362665,	39.99684906],
                [3.325952,	-186.7142944,	3.99222827,	    39.96478271],
                [3.465216,	-186.7851563,	3.991145134,	39.96478271],
                [3.620864,	-186.8465576,	3.990140915,	39.96478271],
                [3.735552,	-186.8771973,	3.989127159,	39.96478271],
                [3.866624,	-186.940918,	3.988111734,	39.96478271],
                [3.997696,	-186.9778442,	3.987061501,	39.96478271],
                [4.161536,	-187.0187988,	3.985952377,	40.02095795],
                [4.3008,	-187.1199951,	3.984705925,	40.02095795],
                [4.46464,	-187.1402588,	3.983659267,	40.02095795],
                [4.612096,	-187.1956177,	3.982604742,	40.02095795],
                [4.759552,	-187.2490234,	3.981431246,	40.10401917],
                [4.939776,	-187.2973633,	3.980323792,	40.10401917],
                [5.095424,	-187.3647461,	3.979264975,	40.10401917],
                [5.24288,	-187.4018555,	3.978199244,	40.10401917],
                [5.431296,	-187.4415894,	3.977194071,	40.10401917],
                [5.603328,	-187.510376,	3.976042271,	40.02743149],
                [5.767168,	-187.5344238,	3.974952221,	40.02743149]])


para_set = pb.ParameterValues(para)

curr_interp_function = pb.Interpolant(dc[:, 0], -dc[:, 1], pb.t)
T_interp_function = pb.Interpolant(dc[:, 0], dc[:, 3]+273.15, pb.t)

para_set["Current function [A]"] = curr_interp_function
para_set["Initial temperature [K]"] = dc[0, 3] + 273.15
para_set["Ambient temperature [K]"] = dc[0, 3] + 273.15
# para_set["Ambient temperature [K]"] = T_interp_function

options = {"cell geometry": "arbitrary", "operating mode": "current"}
mymodel = pb.lithium_ion.DFN(options = options)
var = pb.standard_spatial_vars
var_pts = {var.x_n: 50, var.x_s: 15, var.x_p: 50, var.r_n: 20, var.r_p: 20}

solver_opt = pb.CasadiSolver(mode = 'safe', dt_max = 1e-6)

sim = pb.Simulation(mymodel, parameter_values = para_set, solver = solver_opt, var_pts = var_pts)
t_eval = np.arange(0, 5, 0.1)
sol = sim.solve(t_eval)

time = sol["Time [s]"].data
volt = sol["Terminal voltage [V]"].data

plt.plot(time, volt)

swhan22 Nov 28, 2023
Author

it appears that with v23.5, I am able to set ["Ambient temperature [K]"] as pybamm.Interpolant without getting any error. Some reason, v23.9 gives me the broadcast error. I'll stick with v23.5 for now. Thank you.

rtimms Nov 28, 2023
Collaborator

Ambient temperature can now depend on y,z,t, so you'll need to wrap your interpolant function. Somthing like

def T_function(y,z,t):
    return T_interp_function(t)

Answer selected by swhan22

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Uh oh!

setting ["Ambient temperature [K]"] as pybamm.Interpolant and getting broadcast error #3568

Uh oh!

{{title}}

Uh oh!

Replies: 1 comment 4 replies

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{editor}}'s edit

{{editor}}'s edit

Uh oh!

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{title}}

Uh oh!

Select a reply

Uh oh!

Uh oh!

setting ["Ambient temperature [K]"] as pybamm.Interpolant and getting broadcast error #3568

Uh oh!

swhan22 Nov 27, 2023

Replies: 1 comment · 4 replies

Uh oh!

valentinsulzer Nov 27, 2023 Maintainer

Uh oh!

Uh oh!

swhan22 Nov 27, 2023 Author

Uh oh!

swhan22 Nov 28, 2023 Author

Uh oh!

rtimms Nov 28, 2023 Collaborator

Uh oh!

rtimms Nov 28, 2023 Collaborator

swhan22
Nov 27, 2023

Replies: 1 comment 4 replies

valentinsulzer
Nov 27, 2023
Maintainer

swhan22 Nov 27, 2023
Author

swhan22 Nov 28, 2023
Author

rtimms Nov 28, 2023
Collaborator

rtimms Nov 28, 2023
Collaborator