Source code for mtnlion.models.isothermal

"""
Base Newman isothermal model
"""
import dolfin as fem
import ufl  # type: ignore

import mtnlion.formulas.approximation
import mtnlion.formulas.dfn
from mtnlion.formula import Formula
from mtnlion.formulas.approximation import LagrangeMultiplier
from mtnlion.model import Model
from mtnlion.variable import Variable

# pylint: disable=invalid-name


class Isothermal(Model):
    """
    The basic DFN lithium-ion model with no thermal considerations and a 1D approximation of the solid concentration.
    """

    def __init__(self, num_functions):
        """
        Defines the basic Newman Isothermal model.

        :param num_functions: Number of functions to approximate `self.SolidConcentration`
        """
        super(Isothermal, self).__init__(minimum_rothes_order=1)

        self.Ns = num_functions
        self.legendre = mtnlion.formulas.approximation.Legendre(self.Ns)

        self.variables = [
            Variable("j", ["anode", "cathode"]),
            Variable("phi_s", ["anode", "cathode"]),
            Variable("phi_e", ["anode", "separator", "cathode"]),
            Variable("c_s", ["anode", "cathode"], self.Ns),
            Variable("c_e", ["anode", "separator", "cathode"]),
            Variable("lm_phis_gnd", ["anode_cc"]),
            Variable("lm_phie", ["anode-separator", "separator-cathode"]),
            Variable("lm_ce", ["anode-separator", "separator-cathode"]),
        ]

        self.formulas = [
            self.IntercalationFlux(),
            self.SolidPotential(),
            self.ElectrolytePotential(),
            self.SolidConcentration(self.legendre),
            self.ElectrolyteConcentration(),
            self.SolidConcentrationBoundary(),
            self.SolidPotentialNeumann(),
            self.SolidConcentrationNeumann(),
            self.ExchangeCurrentDensity(),
            self.Overpotential(),
            LagrangeMultiplier(["anode_cc"], "lm_phis_gnd", "phi_s"),
            LagrangeMultiplier(["anode-separator", "separator-cathode"], "lm_phie", "phi_e"),
            LagrangeMultiplier(["anode-separator", "separator-cathode"], "lm_ce", "c_e"),
            self.StateOfCharge(),
            self.OpenCircuitPotential(),
            self.KappaRef(),
            self.KappaEff(),
            self.KappaDEff(),
        ]

    class IntercalationFlux(Formula):
        """
        Describes how the electrical current on an electrode depends on the electrode potential.
        """

        def __init__(self):
            super(Isothermal.IntercalationFlux, self).__init__(domains=["anode", "cathode"])
            self.Variables = self.typedef("Variables", "j")
            self.Parameters = self.typedef("Parameters", "alpha, F, R, T, i0, eta")

        def form(self, arguments, domain):
            (j,) = arguments.variables
            alpha, F, R, T, i0, eta = arguments.parameters

            return (
                j - i0 * (fem.exp((1 - alpha) * F * eta / (R * T)) - fem.exp(-alpha * F * eta / (R * T)))
            ) * j.test()

    class SolidPotential(Formula):
        """
        Charge conservation in the solid.
        """

        def __init__(self):
            super(Isothermal.SolidPotential, self).__init__(domains=["anode", "cathode"])
            self.Variables = self.typedef("Variables", "phi_s, j")
            self.Parameters = self.typedef("Parameters", "a_s, F, sigma_eff, L")

        def form(self, arguments, domain):
            phi_s, j = arguments.variables
            a_s, F, sigma_eff, L = arguments.parameters

            lhs = -sigma_eff / L * fem.dot(fem.grad(phi_s.trial(0)), fem.grad(phi_s.test()))
            rhs = L * a_s * F * j * phi_s.test()

            return rhs - lhs

    class ElectrolytePotential(Formula):
        """
        Charge conservation in the electrolyte.
        """

        def __init__(self):
            super(Isothermal.ElectrolytePotential, self).__init__(domains=["anode", "separator", "cathode"])
            self.Variables = self.typedef("Variables", "phi_e, c_e, j")
            self.Parameters = self.typedef("Parameters", "L, a_s, F, kappa_eff, kappa_Deff")
            self.Formulas = self.typedef("Formulas", "")

        def form(self, arguments, domain):
            phie, ce, j = arguments.variables
            L, a_s, F, kappa_eff, kappa_Deff = arguments.parameters

            if domain == "separator":
                j = 0

            lhs = kappa_eff / L * fem.dot(fem.grad(phie.trial(0)), fem.grad(phie.test()))  # all domains
            rhs1 = -kappa_Deff / L * fem.dot(fem.grad(fem.ln(ce.trial(0))), fem.grad(phie.test()))  # all domains
            rhs2 = L * a_s * F * j * phie.test()  # electrodes

            return lhs - rhs1 - rhs2

    class SolidConcentration(Formula):
        """
        Concentration of lithium in the solid, 1D approximation using Legendre polynomials.
        """

        def __init__(self, legendre):
            super(Isothermal.SolidConcentration, self).__init__(domains=["anode", "cathode"])
            self.Variables = self.typedef("Variables", "c_s")
            self.TimeDiscretization = self.typedef("TimeDiscretization", "dt_cs")
            self.Parameters = self.typedef("Parameters", "Rs, Ds_ref")
            self.legendre = legendre

        def form(self, arguments, domain):
            (cs,) = arguments.variables
            Rs, Ds_ref = arguments.parameters
            (dt_cs,) = arguments.time_discretization

            Ns = self.legendre.num_functions
            M = self.legendre.M
            K = self.legendre.K

            lhs_terms = [M[m, n] * dt_cs(cs.trial(subfunction=n)) * cs.test(m) for m in range(Ns) for n in range(Ns)]
            rhs_terms = [K[m, n] * cs.trial(0, n) * cs.test(m) for m in range(Ns) for n in range(Ns)]

            lhs = sum(lhs_terms) * Rs
            rhs = sum(rhs_terms) * Ds_ref / Rs

            return rhs + lhs

    class ElectrolyteConcentration(Formula):
        """
        Concentration of lithium in the electrolyte.
        """

        def __init__(self):
            super(Isothermal.ElectrolyteConcentration, self).__init__(domains=["anode", "separator", "cathode"])
            self.Variables = self.typedef("Variables", "c_e, j")
            self.Parameters = self.typedef("Parameters", "a_s, De_eff, t_plus, L, eps_e")
            self.TimeDiscretization = self.typedef("TimeDiscretization", "dt_ce")

        def form(self, arguments, domain):
            ce, j = arguments.variables
            a_s, De_eff, t_plus, L, eps_e = arguments.parameters
            (dt_ce,) = arguments.time_discretization

            if domain == "separator":
                j = 0

            lhs = L * eps_e * ce.test()  # all domains
            rhs1 = -De_eff / L * fem.dot(fem.grad(ce.trial(0)), fem.grad(ce.test()))  # all domains
            rhs2 = L * a_s * (1 - t_plus) * j * ce.test()  # electrodes

            return lhs * dt_ce(ce.trial()) - rhs1 - rhs2

    class SolidConcentrationBoundary(Formula):
        """
        This `Formula` defines the value of the lithium concentration at the surface of the solid particle.
        """

        def __init__(self):
            super(Isothermal.SolidConcentrationBoundary, self).__init__(name="cse", domains=["anode", "cathode"])
            self.Variables = self.typedef("Variables", "c_s")

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

            return sum(cs)

    class SolidPotentialNeumann(Formula):
        """
        Neumann boundary for the solid potential at the anode/cathode current collector boundaries.
        """

        def __init__(self):
            super(Isothermal.SolidPotentialNeumann, self).__init__(domains=["anode_cc", "cathode_cc"])
            self.Variables = self.typedef("Variables", "phi_s")
            self.Parameters = self.typedef("Parameters", "Iapp, Acell")

        def form(self, arguments, domain):
            (phis,) = arguments.variables
            Iapp, Acell = arguments.parameters

            return Iapp / Acell * phis.test()

    class SolidConcentrationNeumann(Formula):
        """
        Nuemann boundary for the solid concentration. This `Formula` doesn't use a boundary domain due to the 1D
        1D approximation.
        """

        def __init__(self):
            super(Isothermal.SolidConcentrationNeumann, self).__init__(domains=["anode", "cathode"])
            self.Variables = self.typedef("Variables", "c_s, j")

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

            return j * sum(cs.test(all_funcs=True))

    class ExchangeCurrentDensity(Formula):
        """
        The exchange current density is the current in the absence of net electrolysis and at zero overpotential.
        """

        def __init__(self):
            super(Isothermal.ExchangeCurrentDensity, self).__init__(name="i0", domains=["anode", "cathode"])
            self.Variables = self.typedef("Variables", "c_e")
            self.Parameters = self.typedef("Parameters", "csmax, ce0, alpha, k_norm_ref, cse")

        def form(self, arguments, domain):
            (ce,) = arguments.variables
            csmax, ce0, alpha, k_norm_ref, cse = arguments.parameters

            return (
                k_norm_ref
                * (abs(((csmax - cse) / csmax)) ** (1 - alpha))
                * (abs(cse / csmax) ** alpha)
                * (abs(ce / ce0) ** (1 - alpha))
            )

    class Overpotential(Formula):
        """
        Voltage difference between a reduction potential and the potential of the redox event.
        """

        def __init__(self):
            super(Isothermal.Overpotential, self).__init__(name="eta", domains=["anode", "cathode"])
            self.Variables = self.typedef("Variables", "phi_s, phi_e")
            self.Parameters = self.typedef("Parameters", "Uocp")

        def form(self, arguments, domain):
            phis, phie = arguments.variables
            (Uocp,) = arguments.parameters

            return phis - phie - Uocp

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

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

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

            return cse / csmax

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

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

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

            if domain == "anode":
                return -0.16 + 10.0 * fem.exp(-2000.0 * soc) + 1.32 * fem.exp(-3.0 * soc)

            if domain == "cathode":
                return (
                    -0.0275479 * (0.998432 - soc) ** (-0.492465)
                    + 0.0565661 * ufl.tanh(8.60942 - 14.5546 * soc)
                    + 4.250661588169
                    - 0.157123 * fem.exp(-0.04738 * soc ** 8)
                    + 171.498408536192 * fem.exp(-40 * soc)
                )

            return None

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

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

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

            return -1.6018e-14 * c_e ** 4 + 1.5094e-10 * c_e ** 3 - 4.7212e-7 * c_e ** 2 + 0.0005007 * c_e + 0.041253

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

        def __init__(self):
            super(Isothermal.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(Isothermal.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