Source code for mtnlion.formulas.dfn

"""
A collection of formulas useful for the Doyl-Fuller-Newman cell model
"""
import dolfin as fem
import ufl  # type: ignore
from ufl import algebra  # type: ignore
from ufl.core.expr import Expr  # type: ignore

from mtnlion.domain import Domain
from mtnlion.formula import Arguments, Formula

# pylint: disable=W0123, C0103

SAFE_DICT = {"abs": algebra.Abs}
FENICS_SAFE_LIST = [
    "max_value",
    "min_value",
    "sign",
    "sqrt",
    "exp",
    "ln",
    "erf",
    "cos",
    "sin",
    "tan",
    "acos",
    "asin",
    "atan",
    "atan_2",
    "cosh",
    "sinh",
    "tanh",
    "bessel_J",
    "bessel_Y",
    "bessel_I",
    "bessel_K",
]

SAFE_DICT.update({(k, getattr(ufl, k)) for k in FENICS_SAFE_LIST})


def eval_form(formula: Formula, *parameters) -> Domain[str, Expr]:
    """
    Evaluate a formula with the provided parameters. Note: The order of the parameters must match the definition of the
    formula.

    :param formula: Formula to evaluate
    :param parameters: Parameters to use for evaluation
    :return: FFL expression
    """
    output = {}
    for domain in formula.domains:
        arguments = Arguments(parameters=(parameter[domain] for parameter in parameters))
        output[domain] = formula.formulate(arguments, domain)

    return Domain(output)


class SOC(Formula):
    """State of Charge (SOC) formula."""

    def __init__(self):
        super(SOC, self).__init__(name="soc", domains=["anode", "cathode"])
        self.Parameters = self.typedef("Parameters", "cse, csmax")

    def form(self, arguments, domain):
        (cse, csmax,) = arguments.parameters

        return cse / csmax


class Uocp(Formula):
    """Open-circuit potential formula."""

    def __init__(self, uocp_str):
        super(Uocp, self).__init__(name="Uocp", domains=["anode", "cathode"])
        self.Parameters = self.typedef("Parameters", "soc")
        self.uocp_str = uocp_str

    def form(self, arguments, domain):
        """Evaluate the open-circuit potential equation."""
        (soc,) = arguments.parameters

        return eval(str(self.uocp_str[domain]), {"__builtins__": None}, {"soc": soc, **SAFE_DICT})


class KappaRef(Formula):
    """Bulk conductivity of the homogeneous materials."""

    def __init__(self, kappa_ref_str):
        super(KappaRef, self).__init__(name="kappa_ref", domains=["anode", "separator", "cathode"])
        self.Variables = self.typedef("Variables", "c_e")
        self.kappa_ref_str = kappa_ref_str

    def form(self, arguments, domain):
        (c_e,) = arguments.variables

        return eval(str(self.kappa_ref_str), {"__builtins__": None}, {"x": c_e, **SAFE_DICT})


class KappaEff(Formula):
    """Effective conductivity of the electrolyte."""

    def __init__(self):
        super(KappaEff, self).__init__(name="kappa_eff", domains=["anode", "separator", "cathode"])
        self.Parameters = self.typedef("Parameters", "eps_e, brug_kappa, kappa_ref")

    def form(self, arguments, domain):
        eps_e, brug_kappa, kappa_ref = arguments.parameters

        return kappa_ref * eps_e ** brug_kappa


class KappaDEff(Formula):
    """kappa_d effective"""

    def __init__(self):
        super(KappaDEff, self).__init__(name="kappa_Deff", domains=["anode", "separator", "cathode"])
        self.Parameters = self.typedef("Parameters", "eps_e, kappa_D, kappa_ref")

    def form(self, arguments, domain):
        eps_e, kappa_d, kappa_ref = arguments.parameters

        return kappa_d * kappa_ref * eps_e


class FilmThickness(Formula):
    """
    Thickness of the film that builds around the surface of the particles.
    """

    def __init__(self):
        super(FilmThickness, self).__init__(domains=["anode"])
        self.Variables = self.typedef("Variables", "delta_film, j_s")
        self.Parameters = self.typedef("Parameters", "Mp, rho_p")
        self.TimeDiscretization = self.typedef("TimeDiscretization", "dt_delta_film")

    def form(self, arguments, domain):
        delta_film, js = arguments.variables
        Mp, rho_p = arguments.parameters
        (dt,) = arguments.time_discretization

        lhs = dt(delta_film.trial()) * delta_film.test()
        rhs = -Mp / rho_p * js * delta_film.test()
        return lhs - rhs


class FilmResistance(Formula):
    """
    Resistance of the film that builds around the surface of the particles.
    """

    def __init__(self):
        super(FilmResistance, self).__init__(name="Rfilm", domains=["anode"])
        self.Variables = self.typedef("Variables", "delta_film")
        self.Parameters = self.typedef("Parameters", "Rsei, kappa_p")

    def form(self, arguments, domain):
        (delta_film,) = arguments.variables
        Rsei, kappa_p = arguments.parameters

        return Rsei + delta_film / kappa_p


class CapacityLoss(Formula):
    """
    Capacity loss due to the side reactions.
    """

    def __init__(self):
        super(CapacityLoss, self).__init__(domains=["anode"])
        self.Variables = self.typedef("Variables", "Q, j_s")
        self.Parameters = self.typedef("Parameters", "a_s, Acell, F, dx")
        self.TimeDiscretization = self.typedef("TimeDiscretization", "dt_Q")

    def form(self, arguments, domain):
        Q, js = arguments.variables
        a_s, Acell, F, dx = arguments.parameters
        (dt,) = arguments.time_discretization

        lhs = dt(Q.trial()) * Q.test()
        rhs = fem.assemble(a_s * Acell * F * js * dx) * Q.test()
        return lhs - rhs