Create an exchange current density function depends on user-defined parameter.

[samwong92]

I am trying to set the exchange current density to a tanh function of stoichiometry, and the constants (Fitting param1~3) needed to be set with param.update() . (Because I am doing auto fitting with bayes-opt. )
I followed https://github.com/pybamm-team/PyBaMM/blob/develop/examples/notebooks/Creating%20Models/5-half-cell-model.ipynb and some other examples and tried with the code below.
When I plot the exchange current density of the solution, its value is far from I expected, and seems depending on the wrong variables. Please points out my mistake and advice.

import os
import matplotlib.pyplot as plt
import numpy as np
import pybamm
import pandas as pd
import timeit
from matplotlib import style
from pybamm import tanh

model = pybamm.lithium_ion.DFN({"open-circuit potential": ("current sigmoid", "current sigmoid")})
c_surf = pybamm.Variable("R-averaged negative particle concentration [mol.m-3]")
c_max = pybamm.Parameter("Maximum concentration in negative electrode [mol.m-3]")
p1 = pybamm.Parameter("Fitting param1")
p2 = pybamm.Parameter("Fitting param2")
p3 = pybamm.Parameter("Fitting param3")
inputs = {"R-averaged negative particle concentration [mol.m-3]": c_surf,
          "Fitting param1": p1,
          "Fitting param2": p2,
          "Fitting param3": p3
          }
j0 =  pybamm.FunctionParameter("negative electrode exchange-current density [A.m-2]", inputs)
def exchange_current_density_func(c_surf,p1,p2,p3):
    return p1 * tanh(p2*(p3-(c_surf/c_max))) + p1 

param = pybamm.ParameterValues("NewParam2304")
geometry = model.default_geometry
param.update({"Current function [A]": "[input]",
              "Positive electrode diffusivity [m2.s-1]": "[input]",
              "Negative electrode diffusivity [m2.s-1]": "[input]",
              "Positive electrode exchange-current density [A.m-2]": "[input]",
              "Negative electrode exchange-current density [A.m-2]": exchange_current_density_func,
              "Fitting param1": "[input]",
              "Fitting param2": "[input]",
              "Fitting param3": "[input]"
              }, check_already_exists=False)
              
param.process_geometry(geometry)
param.process_model(model)
var = pybamm.standard_spatial_vars
var_pts = {
    "x_n": 40,  # negative electrode
    "x_s": 15,  # separator
    "x_p": 70,  # positive electrode
    "r_n": 20,  # negative particle
    "r_p": 40,  # positive particle
}
mesh = pybamm.Mesh(geometry,model.default_submesh_types, var_pts)
disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
disc.process_model(model)
safe_solver = pybamm.CasadiSolver(atol=1e-3, rtol=1e-3, mode="safe")

solution = safe_solver.solve(
    model, t_eval,
      inputs={"Current function [A]": 1.,
              "Negative electrode diffusivity [m2.s-1]": 9.12e-12,
              "Positive electrode diffusivity [m2.s-1]": 3.04e-12,
              "Positive electrode exchange-current density [A.m-2]": 0.0123,
              "Negative electrode exchange-current density [A.m-2]": exchange_current_density_func,
              "Fitting param1": 0.15,
              "Fitting param2": 10,
              "Fitting param3": 0.1
              }
  )
output_variables = ["X-averaged negative electrode exchange current density [A.m-2]","X-averaged positive electrode exchange current density [A.m-2]","Negative particle surface concentration [mol.m-3]","R-averaged negative particle concentration [mol.m-3]"]
solution.plot(output_variables)

[rtimms]

The model expects the function for the exchange current density to depend on c_e, c_s_surf,c_s_max and T. You can define your function and update the parameter like this

p1 = pybamm.InputParameter("p1")
p2 = pybamm.InputParameter("p2")
p3 = pybamm.InputParameter("p3")

def j0(c_e, c_s_surf, c_s_max, T):
    return p1 * tanh(p2*(p3-(c_surf/c_max))) + p1


param.update({"Negative electrode exchange-current density [A.m-2]" j0})

Then when you solve give values for "p1" etc. in inputs.

[samwong92]

Thank you, that solved my problem. I have a further question hope you can answer me.
Could I modify the electrode diffusivity to be stoichiometry-dependent using the same method?
For example,

def D_p_sto(sto):
    return p1 * tanh(p2*(p3-(sto))) + p1

param.update({"Positive electrode diffusivity [m2.s-1]": D_p_sto})

[samwong92]

Let me reply myselves.
Depend on what default model I am using, e.g. basic_dfn.py,
:169 ~ 170

        N_s_n = -param.n.prim.D(c_s_n, T) * pybamm.grad(c_s_n)
        N_s_p = -param.p.prim.D(c_s_p, T) * pybamm.grad(c_s_p)

The model pass (c_s, T) to the diffusivity function.
So if I want to modify its dependence to e.g. sto, I need to change those lines too.

Okay, I find that parameters arguments are defined in parameter/____.py , and diffusivity is depends on sto on default, so, let just say I want to add stoichimetry-dependent bruggeman constant, b_s or b_e. I believed change needed to be made here. I wonder is there a way to make changes without messing up the original code? incase of version up etc. maybe a submodel ? building a submodel from zero seems difficult. Any guide could help?

[rtimms]

PyBaMM expects diffusivity to be given as function of sto and T, so do

def D_p_sto(sto, T):
    return p1 * tanh(p2*(p3-(sto))) + p1

param.update({"Positive electrode diffusivity [m2.s-1]": D_p_sto})

As you say, you can find all the arguments in parameters/*.py (e.g. see [here](https://github.com/pybamm-team/PyBaMM/blob/develop/pybamm/parameters/lithium_ion_parameters.py#L521 for diffusivity). If you need to change these you will need to change the source code to call the parameters with the new arguments.

[samwong92]

Thank you @rtimms.
I have a problem when modifying stoichimetry-dependent Bruggeman constant for electrode, b_e, in pybamm/parameters/lithium_ion_parameters.py.
I mimic D(self, c_s, T) and defined a new function to set b_e as a functionParameter,
and changed basic_dfn.py as below

        tor = pybamm.concatenation(
          # eps_n**param.n.b_e, eps_s**param.s.b_e, eps_p**param.p.b_e
            eps_n**param.n.b_e(c_s_n), eps_s**param.s.b_e, eps_p**param.p.b_e(c_s_p) 
            )

However, symbol errors occur when I tried updating b_e as a function.
I believe it is because the defined parameterFunction not being used. Then, If I delete self.geo.b_e in _set_parameters(self), " "NotImplementedError",
"message": "BinaryOperator not implemented for symbols of type <class 'pybamm.expression_tree.parameter.FunctionParameter'> and <class 'method'>"
occurs.
Since no parameterFunction is defined inside _set_parameters(self), I believe this is not the solution.
Any advice?
Below is my b_c

def b_e(self, c_s):
     Domain = self.domain.capitalize()
     sto = c_s / self.c_max
     tol = pybamm.settings.tolerances["D__c_s"]
     sto = pybamm.maximum(pybamm.minimum(sto, 1 - tol), tol)
     inputs = {
         f"{Domain} particle stoichiometry": sto
         # f"{self.phase_prefactor}{Domain} particle stoichiometry": sto,
         # "Temperature [K]": T,
     }
     return pybamm.FunctionParameter(
         f"{Domain} electrode Bruggeman coefficient (electrolyte)",
         inputs,
     )