Adding new variable to Positive OCP Voltage parameter function

[pysam305]

Hi, I am trying to introduce a concentration term in my Positive OCP function.
The code is pasted below. Sorry for the raw format

import pybamm
import os
pybamm.set_logging_level("NOTICE")
model=pybamm.lithium_ion.DFN(options={"SEI":"solvent-diffusion limited","particle mechanics":"swelling only"})
params=pybamm.ParameterValues("OKane2022")
def nmc_LGM50_ocp_Chen2020(sto,c_s_max):
    """
    LG M50 NMC open-circuit potential as a function of stochiometry, fit taken
    from [1].
    References
    ----------
    .. [1] Chang-Hui Chen, Ferran Brosa Planella, Kieran O’Regan, Dominika Gastol, W.
    Dhammika Widanage, and Emma Kendrick. "Development of Experimental Techniques for
    Parameterization of Multi-scale Lithium-ion Battery Models." Journal of the
    Electrochemical Society 167 (2020): 080534.
    Parameters
    ----------
    sto: :class:`pybamm.Symbol`
        Electrode stochiometry
    Returns
    -------
    :class:`pybamm.Symbol`
        Open-circuit potential
    """

    u_eq = (
        -0.8090 * sto
        + 4.4875
        - 0.0428 * pybamm.tanh(18.5138 * (sto - 0.5542))
        - 17.7326 * pybamm.tanh(15.7890 * (sto - 0.3117))
        + 17.5842 * pybamm.tanh(15.9308 * (sto - 0.3120))
    )
    E=params["Positive electrode Young's modulus [Pa]"]
    omega = params["Positive electrode partial molar volume [m3.mol-1]"]
    poisson=params["Positive electrode Poisson's ratio"]
    gamma_stress=(E*omega*omega)/(3*(1-poisson))*(c_s_max/pybamm.constants.F)

    return u_eq+gamma_stress
params["Positive electrode OCP [V]"]=nmc_LGM50_ocp_Chen2020(model.variables["Positive particle stoichiometry"],model.variables["Maximum positive particle concentration [mol.m-3]"])
exp=pybamm.Experiment([("Rest for 3600 seconds","Charge at 1.6252 A until 4.2 V","Hold at 4.2 V until 0.05 A","Rest for 1 hour","Discharge at 1 C until 2.5V","Rest for 1 hour")]*1)
solver = pybamm.CasadiSolver(mode="safe",dt_max=1e-150,max_step_decrease_count=1e-1)
sim=pybamm.Simulation(model=model,experiment=exp, parameter_values=params,solver=solver)
sim.solve()

The error being shown here is

"Parameter provided for 'Positive electrode OCP [V]' is of the wrong type (should either be scalar-like or callable)"

Could anyone please help as to how I can get over this. Also I don't understand how without mentioning what sto is how pybamm runs the code, like, if I just write
"params["Positive electrode OCP [V]"]=nmc_LGM50_ocp_Chen2020" and not include c_s_max the simulation runs just fine but the moment I include something like "T" in the nmc function as argument it say
"TypeError: nmc_LGM50_ocp_Chen2020() missing 1 required positional argument: 'T'".
How does it run fine without sto but not with addition arguments.
Could anyone please solve this dilemma. Thanks.

[brosaplanella]

OCP should only be a function of the stoichiometry, so if you add any other arguments it will throw an error. If you need c_s_max you can get it through the parameters, similarly as with omega and poisson.

[pysam305]

Hi. Thanks for the reply.
I have followed your suggestion and used c_s_max from the parameter list. However, when specifying it the following error came.
def nmc_LGM50_ocp_Chen2020(sto):
"""
LG M50 NMC open-circuit potential as a function of stochiometry, fit taken
from [1].
References
----------
.. [1] Chang-Hui Chen, Ferran Brosa Planella, Kieran O’Regan, Dominika Gastol, W.
Dhammika Widanage, and Emma Kendrick. "Development of Experimental Techniques for
Parameterization of Multi-scale Lithium-ion Battery Models." Journal of the
Electrochemical Society 167 (2020): 080534.
Parameters
----------
sto: :class:pybamm.Symbol
Electrode stochiometry
Returns
-------
:class:pybamm.Symbol
Open-circuit potential
"""

u_eq = (
    -0.8090 * sto
    + 4.4875
    - 0.0428 * pybamm.tanh(18.5138 * (sto - 0.5542))
    - 17.7326 * pybamm.tanh(15.7890 * (sto - 0.3117))
    + 17.5842 * pybamm.tanh(15.9308 * (sto - 0.3120))
)
E=params["Positive electrode Young's modulus [Pa]"]
omega = params["Positive electrode partial molar volume [m3.mol-1]"]
poisson=params["Positive electrode Poisson's ratio"]
c_s_max=params["Maximum concentration in positive electrode [mol.m-3]"]
gamma_stress=(E*omega*omega)/(3*(1-poisson))*(c_s_max/pybamm.constants.F)

return u_eq+gamma_stress

params["Positive electrode OCP [V]"]=nmc_LGM50_ocp_Chen2020
exp=pybamm.Experiment([("Rest for 3600 seconds","Charge at 1.6252 A until 4.2 V","Hold at 4.2 V until 0.05 A","Rest for 1 hour","Discharge at 1 C until 2.5V","Rest for 1 hour")]*1)
sim=pybamm.Simulation(model=model,experiment=exp, parameter_values=params,solver=solver)
sim.solve()

The error being shown is:
Events ['Maximum voltage [V]'] are non-positive at initial conditions
Now I tried to change the initial SOC but whatever I do it says reduce dt_max.
sim=pybamm.Simulation(model=model,experiment=exp, parameter_values=params,solver=solver)
sim.solve(initial_soc=0.1)

The above code is for dt_max=1e-200,max_step_decrease_count=1e-10, and the error being shown is

Maximum number of decreased steps occurred at t=0.0 (final SolverError: 'Error in Function::call for 'F' [IdasInterface] at .../casadi/core/function.cpp:1401:
Error in Function::call for 'F' [IdasInterface] at .../casadi/core/function.cpp:330:
.../casadi/interfaces/sundials/idas_interface.cpp:596: IDACalcIC returned "IDA_LINESEARCH_FAIL". Consult IDAS documentation.'). For a full solution try reducing dt_max (currently, dt_max=60.6) and/or reducing the size of the time steps or period of the experiment. Set return_solution_if_failed_early=True to return the solution object up to the point where failure occured.

Could anyone please tell me how to get over this roadblock