#!/usr/bin/env python
# -*- coding: utf-8 -*-

                Copyright (C) Raymond B. Smith, 2016
This program is free software; you can redistribute it and/or modify it under the
terms of the GNU General Public License version 3 as published by the Free Software
Foundation. This program is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
PARTICULAR PURPOSE. See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License along with the
DAE Tools software; if not, see <http://www.gnu.org/licenses/>.
from daetools.pyDAE import *
import numpy as np

from pyUnits import m, s, K, mol, J, A, V, S

__doc__ = """
Model of a lithium-ion battery based on porous electrode theory as developed
by John Newman and coworkers. In particular, the equations here are based on a summary
of the methodology by Karen E. Thomas, John Newman, and Robert M. Darling,

Thomas K., Newman J., Darling R. (2002). Mathematical Modeling of Lithium Batteries
in Advances in Lithium-ion Batteries. Springer US. 345-392.
`doi:10.1007/0-306-47508-1_13 <http://dx.doi.org/10.1007/0-306-47508-1_13>`_

A few simplifications have been made rather than implementing the more complete model described there.
For example, the following assumptions have (currently) been made:

- two porous electrodes are used rather than providing the option for a "half cell" in which
  one electrode is lithium foil.
- conductivity in the electron-conducting phase is infinite
- constant exchange current density in Butler-Volmer reaction expression
- no electrolyte convection
- constant and uniform solvent concentration (ions vary according to concentrated solution theory)
- monodisperse particles in electrode
- no volume occupied by binder, filler, etc. in the electrode


# Define some variable types
conc_t = daeVariableType(
    name="conc_t", units=mol/(m**3), lowerBound=0,
    upperBound=1e20, initialGuess=1.00, absTolerance=1e-6)
elec_pot_t = daeVariableType(
    name="elec_pot_t", units=V, lowerBound=-1e20,
    upperBound=1e20, initialGuess=0, absTolerance=1e-5)
current_dens_t = daeVariableType(
    name="current_dens_t", units=A/m**2, lowerBound=-1e20,
    upperBound=1e20, initialGuess=0, absTolerance=1e-5)
rxn_t = daeVariableType(
    name="rxn_t", units=mol/(m**2 * s), lowerBound=-1e20,
    upperBound=1e20, initialGuess=0, absTolerance=1e-5)

def kappa(c):
    """Return the conductivity of the electrolyte in S/m as a function of concentration in M."""
    out = 0.1  # S/m
    return out * Constant(1 * S/m)

def D(c):
    """Return electrolyte diffusivity (in m^2/s) as a function of concentration in M."""
    out = 1e-10  # m**2/s
    return out * Constant(1 * m**2/s)

def thermodynamic_factor(c):
    """Return the electrolyte thermodynamic factor as a function of concentration in M."""
    out = 1
    return out

def t_p(c):
    """Return the electrolyte cation transference number as a function of concentration in M."""
    out = 0.3 * c/c
    return out

def Ds_n(y):
    """Return diffusivity (in m^2/s) as a function of solid filling fraction, y."""
    out = 5e-9 * y/y  # m**2/s
    return out * Constant(1 * m**2/s)

def Ds_p(y):
    """Return diffusivity (in m^2/s) as a function of solid filling fraction, y."""
    out = 1e-9 * y/y  # m**2/s
    return out * Constant(1 * m**2/s)

def U_n(y):
    """Return the equilibrium potential (V vs Li) of the negative electrode active material
    as a function of solid filling fraction, y.
    material = "coke"
    if material == "coke":
        # Carbon (coke) -- Fuller, Doyle, Newman, J. Electrochem. Soc., 1994
        out = -0.132 + 1.42*np.exp(-2.52*(0.5*y))
    elif material == "Li metal":
        # Lithium metal
        out = 0.
    units = Constant(1 * V) if isinstance(y, adouble) else V
    return out * units

def U_p(y):
    """Return the equilibrium potential (V vs Li) of the positive electrode active material
    as a function of solid filling fraction, y.
    material = "Mn2O4"
    if material == "Mn2O4":
        # Mn2O4 -- Fuller, Doyle, Newman, J. Electrochem. Soc., 1994
        out = (4.06279 + 0.0677504*np.tanh(-21.8502*y + 12.8268) -
               0.105734*(1/((1.00167 - y)**(0.379571)) - 1.575994) -
               0.045*np.exp(-71.69*y**8) +
               0.01*np.exp(-200*(y - 0.19)))
    elif material == "Li metal":
        # Lithium metal
        out = 0.
    units = Constant(1 * V) if isinstance(y, adouble) else V
    return out * units

class ModParticle(daeModel):
    def __init__(self, Name, pindx1, pindx2, c2, y_avg, phi2, phi1, Ds, U, Parent=None, Description=""):
        daeModel.__init__(self, Name, Parent, Description)
        self.Ds = Ds
        self.U = U

        # Domain where variables are distributed
        self.r = daeDomain("r", self, m, "radial domain in particle")

        # Variables
        self.c = daeVariable("c", conc_t, self, "Concentration in the solid")
        self.j_p = daeVariable("j_p", rxn_t, self, "Rate of reaction into the solid")

        # Parameter
        self.j_0 = daeParameter("j_0", mol/(m**2 * s), self, "Exchange current density / F")
        self.alpha = daeParameter("alpha", unit(), self, "Reaction symmetry factor")
        self.c_ref = daeParameter("c_ref", mol/m**3, self, "Max conc of species in the solid")
        self.V_thermal = daeParameter("V_thermal", V, self, "Thermal voltage")
        self.R = daeParameter("R", m, self, "Radius of particle")

        self.pindx1 = pindx1
        self.pindx2 = pindx2
        self.phi2 = phi2
        self.c2 = c2
        self.phi1 = phi1
        self.y_avg = y_avg

    def DeclareEquations(self):

        # Mass conservation in the solid particles governed by (possibly non-linear) Ficks Law diffusion
        # Thomas et al., Eq 17
        eq = self.CreateEquation("mass_cons")
        r = eq.DistributeOnDomain(self.r, eOpenOpen)
        c = self.c(r)
        dc = d(c, self.r, eCFDM)
        D = self.Ds(c/self.c_ref())
        eq.Residual = dt(c) - 1/r()**2*d(r()**2*D*dc, self.r, eCFDM)

        # Symmetry at the center from particles with spherical geometry and symmetry
        # Thomas et al., Eq 18
        eq = self.CreateEquation("CenterSymmetry", "dc/dr = 0 at r=0")
        r = eq.DistributeOnDomain(self.r, eLowerBound)
        c = self.c(r)
        eq.Residual = d(c, self.r, eCFDM)

        # Flux at the particle surface given by the electrochemical reaction rate of (di)intercalation
        # Thomas et al., Eq 18
        eq = self.CreateEquation("SurfaceGradient", "D_s*dc/dr = j_+ at r=R_p")
        r = eq.DistributeOnDomain(self.r, eUpperBound)
        c = self.c(r)
        eq.Residual = self.Ds(c/self.c_ref()) * d(c, self.r, eCFDM) - self.j_p()

        # The rate of electrochemical reaction calculated via the Butler-Volmer equation
        # Here, we use a constant exchange current density, but other forms depending on
        # solid and electrolyte concentrations are commonly used.
        # Thomas et al., Eq 19 and 27
        c_surf = self.c(self.r.NumberOfPoints - 1)
        eta = self.phi1(self.pindx1) - self.phi2(self.pindx2) - self.U(c_surf/self.c_ref())
        eta_ndim = eta / self.V_thermal()
        # At time t=0, to aid in initialization, we use a linearized form of the reaction equation
        self.IF(Time() == Constant(0*s))
        eq = self.CreateEquation("SurfaceRxn", "Reaction rate")
        eq.Residual = self.j_p() + self.j_0() * eta_ndim
        # For the rest of the simulation, we use the full Butler-Volmer equation
        eq = self.CreateEquation("SurfaceRxn", "Reaction rate")
        eq.Residual = self.j_p() - self.j_0() * (np.exp(-self.alpha()*eta_ndim) - np.exp((1 - self.alpha())*eta_ndim))

        # For convenience, we also keep track of the average filling fraction in this particle.
        # This is obtained from integrating the conservation equation over the spherical
        # particle and applying the divergence theorem.
        eq = self.CreateEquation("y_avg")
        eq.Residual = self.y_avg.dt(self.pindx2) - 3*self.j_p()/(self.c_ref()*self.R())

class ModCell(daeModel):
    def __init__(self, Name, Parent=None, Description="", process_info=None):
        daeModel.__init__(self, Name, Parent, Description)
        self.process_info = process_info

        # Domains where variables are distributed
        self.x_centers_n = daeDomain("x_centers_n", self, m, "X cell-centers domain in negative electrode")
        self.x_centers_p = daeDomain("x_centers_p", self, m, "X cell-centers domain in positive electrode")
        self.x_centers_full = daeDomain("x_centers_full", self, m, "X cell-centers domain over full cell")
        self.x_faces_full = daeDomain("x_faces_full", self, m, "X cell-faces domain over full cell")

        # Variables
        # Concentration/potential in different regions of electrolyte and electrode
        self.c = daeVariable("c", conc_t, self, "Concentration in the elyte")
        self.phi2 = daeVariable("phi2", elec_pot_t, self, "Electric potential in the elyte")
        self.i2 = daeVariable("i2", current_dens_t, self, "Electrolyte current density")
        self.phi1_n = daeVariable("phi1_n", elec_pot_t, self, "Electric potential in bulk sld, negative")
        self.phi1_p = daeVariable("phi1_p", elec_pot_t, self, "Electric potential in bulk sld, positive")
        self.phiCC_n = daeVariable("phiCC_n", elec_pot_t, self, "phi at negative current collector")
        self.phiCC_p = daeVariable("phiCC_p", elec_pot_t, self, "phi at positive current collector")
        self.y_avg = daeVariable("y_avg", no_t, self, "Average filling fraction in the solid")
        self.V = daeVariable("V", elec_pot_t, self, "Applied voltage")
        self.current = daeVariable("current", current_dens_t, self, "Total current of the cell")

        # Parameters
        self.F = daeParameter("F", A*s/mol, self, "Faraday's constant")
        self.R = daeParameter("R", J/(mol*K), self, "Gas constant")
        self.T = daeParameter("T", K, self, "Temperature")
        self.a_n = daeParameter("a_n", m**(-1), self, "Reacting area per electrode volume, negative electrode")
        self.a_p = daeParameter("a_p", m**(-1), self, "Reacting area per electrode volume, positive electrode")
        self.L_n = daeParameter("L_n", m, self, "Length of negative electrode")
        self.L_s = daeParameter("L_s", m, self, "Length of separator")
        self.L_p = daeParameter("L_p", m, self, "Length of positive electrode")
        self.BruggExp_n = daeParameter("BruggExp_n", unit(), self, "Bruggeman exponent in x_n")
        self.BruggExp_s = daeParameter("BruggExp_s", unit(), self, "Bruggeman exponent in x_s")
        self.BruggExp_p = daeParameter("BruggExp_p", unit(), self, "Bruggeman exponent in x_p")
        self.poros_n = daeParameter("poros_n", unit(), self, "porosity in x_n")
        self.poros_s = daeParameter("poros_s", unit(), self, "porosity in x_s")
        self.poros_p = daeParameter("poros_p", unit(), self, "porosity in x_p")
        self.c_ref = daeParameter("c_ref", mol/m**3, self, "Reference electrolyte concentration")
        self.currset = daeParameter("currset", A/m**2, self, "current per electrode area")
        self.Vset = daeParameter("Vset", V, self, "applied voltage set point")
        self.tau_ramp = daeParameter("tau_ramp", s, self, "Time scale for ramping voltage or current")
        self.xval_cells = daeParameter("xval_cells", m, self, "coordinate of cell centers")
        self.xval_faces = daeParameter("xval_faces", m, self, "coordinate of cell faces")

        # Sub-models
        N_n = self.process_info["N_n"]
        N_s = self.process_info["N_s"]
        N_p = self.process_info["N_p"]
        self.particles_n = np.empty(N_n, dtype=object)
        self.particles_p = np.empty(N_p, dtype=object)
        for indx in range(N_n):
            indx1 = indx2 = indx
            self.particles_n[indx] = ModParticle("particle_n_{}".format(indx), indx1, indx2, self.c,
                                                 self.y_avg, self.phi2, self.phi1_n, Ds_n, U_n, Parent=self)
        for indx in range(N_p):
            indx1 = indx
            indx2 = N_n + N_s + indx
            self.particles_p[indx] = ModParticle("particle_p_{}".format(indx), indx1, indx2, self.c,
                                                 self.y_avg, self.phi2, self.phi1_p, Ds_p, U_p, Parent=self)

    def DeclareEquations(self):

        pinfo = self.process_info
        # Thermal voltage = RT/F = kT/e = approximately 0.026 mV at room temp, 25 C
        V_thm = self.R() * self.T() / self.F()

        # We choose to use (cell centered) finite volume discretization for the electrolyte rather
        # than the built-in finite difference method for a few reasons
        # - It is a mass conservative method, which is important for quasi-neutral electrolyte models
        # - It is more stable at high electrolyte depletion than finite difference methods
        # As a result, we need some information about the domain on which we have discretized our
        # field variables
        # With the finite volume method, we store some information at cell centers
        # (scalar field variables like concentration and electric potential)
        # and store/calculate some information at the faces between cells
        # (fluxes like current density and flux of anions)
        # For more information, see, e.g.,
        #   http://www.ctcms.nist.gov/fipy/documentation/numerical/discret.html
        #   https://en.wikipedia.org/wiki/Finite_volume_method
        # Number of grid cells centers in each of the negative and positive electrodes
        N_n, N_p = self.x_centers_n.NumberOfPoints, self.x_centers_p.NumberOfPoints
        # Number of grid cell centers and faces along the entire electrode
        N_centers, N_faces = self.x_centers_full.NumberOfPoints, self.x_faces_full.NumberOfPoints
        # Number of grid cell centers along the separator
        N_s = N_centers - N_n - N_p
        # Coordinates of the cell centers and cell faces
        center_coords = np.array([self.xval_cells(indx) for indx in range(N_centers)])
        face_coords = np.array([self.xval_faces(indx) for indx in range(N_faces)])
        # Spacing between cell centers, which we will use for finite difference approximations for fluxes
        # We add space for a "ghost point" on each end, which will be used for boundary conditions.
        h_centers = np.hstack((np.diff(center_coords)[0], np.diff(center_coords), np.diff(center_coords)[-1]))
        # Spacing between cell faces.
        h_faces = np.diff(face_coords)

        # For convenience, make numpy arrays of variables at cell centers
        phi2 = np.array([self.phi2(indx) for indx in range(N_centers)])
        c = np.array([self.c(indx) for indx in range(N_centers)])
        dcdt = np.array([self.c.dt(indx) for indx in range(N_centers)])
        a = np.array([self.a_n()]*N_n + [Constant(0 * m**(-1))]*N_s + [self.a_p()]*N_p)
        j_p = np.array([self.particles_n[indx].j_p() for indx in range(N_n)]
                       + [Constant(0 * mol/(m**2 * s))]*N_s
                       + [self.particles_p[indx].j_p() for indx in range(N_p)])
        poros = np.array([self.poros_n()]*N_n + [self.poros_s()]*N_s + [self.poros_p()]*N_p)
        eff_factor_tmp = np.array([self.poros_n() / (self.poros_n()**self.BruggExp_n())]*(N_n+1)
                                  + [self.poros_s() / (self.poros_s()**self.BruggExp_s())]*N_s
                                  + [self.poros_p() / (self.poros_p()**self.BruggExp_p())]*(N_p+1))
        # The eff_factor is a prefactor for the transport in the porous medium compared to transport
        # in a free solution. It is needed at the cell faces because it is used in calculation of fluxes,
        # so we use a harmonic mean to approximate the value at the faces.
        eff_factor = (2*eff_factor_tmp[1:]*eff_factor_tmp[:-1]) / (eff_factor_tmp[1:] + eff_factor_tmp[:-1])

        # Boundary conditions on c and phi2 at current collectors.
        # For concentration, Thomas et al., Eq 15
        # For phi at current collectors, grad(phi) = 0 is required such that i2 = 0 at the current collectors
        # To do these, create "ghost points" on the end of cell-center vectors
        ctmp = np.empty(N_centers + 2, dtype=object)
        ctmp[1:-1] = c
        phi2tmp = np.empty(N_centers + 2, dtype=object)
        phi2tmp[1:-1] = phi2
        # No ionic current passes into the current collectors, which requires
        # grad(c) = grad(phi2) = 0
        # at both current collectors. We apply this by using the ghost points.
        ctmp[0] = ctmp[1]
        ctmp[-1] = ctmp[-2]
        phi2tmp[0] = phi2tmp[1]
        phi2tmp[-1] = phi2tmp[-2]
        # We'll need the value of c at the faces as well. We use a harmonic mean.
        c_faces = (2*ctmp[1:]*ctmp[:-1])/(ctmp[1:] + ctmp[:-1])

        # Approximate the gradients of these field variables at the faces
        dc = np.diff(ctmp) / h_centers
        dlogc = np.diff(np.log(ctmp / self.c_ref())) / h_centers
        dphi2 = np.diff(phi2tmp) / h_centers

        # Effective transport properties are required at faces between cells
        # Thomas et al., below Eq 3
        D_eff = eff_factor * D(c_faces / self.c_ref())
        kappa_eff = eff_factor * kappa(c_faces / self.c_ref())

        # Flux of charge (current density) at faces
        # Thomas et al., Eq 3
        i = -kappa_eff * (dphi2 - 2*V_thm*(1 - t_p(c_faces))*thermodynamic_factor(c_faces)*dlogc)

        # Flux of anions at faces
        # Based on Thomas et al., Eq 8
        # Using Thomas et al. Eq 9 and 10, using z_m = -1 and the paragraph below Eq 12 with Eq 13
        N_m = -D_eff*dc - (1 - t_p(c_faces)) * i / self.F()

        # Store values for the current density
        for indx in range(N_faces):
            eq = self.CreateEquation("i2_{}".format(indx))
            eq.Residual = self.i2(indx) - i[indx]

        # Divergence of fluxes
        di = np.diff(i) / h_faces
        dN_m = np.diff(N_m) / h_faces
        # Electrolyte: mass and charge conservation
        for indx in range(N_centers):
            # Thomas et al., Eq 11
            # Used instead of Eq 12, which is equivalent, noting that c_m = c for the electrolyte
            eq = self.CreateEquation("mass_cons_m_{}".format(indx), "anion mass conservation")
            eq.Residual = poros[indx]*dcdt[indx] + dN_m[indx]
            # Thomas et al., Eq 27 and 28
            eq = self.CreateEquation("charge_cons_{}".format(indx), "charge conservation")
            eq.Residual = -di[indx] - self.F()*a[indx]*j_p[indx]

        # Arbitrary datum for electric potential.
        # Thomas et al., below Eq 3
        # We apply this in the electrolyte at an arbitrary location, the negative current collector
        eq = self.CreateEquation("phi2_datum")
        eq.Residual = self.phiCC_n()

        # Electrode: charge conservation
        phi1_n = np.array([self.phi1_n(indx) for indx in range(N_n)])
        phi1_p = np.array([self.phi1_p(indx) for indx in range(N_p)])
        # We assume infinite conductivity in the electron conducting phase for simplicity
        # negative
        for indx in range(N_n):
            eq = self.CreateEquation("phi1_n_{}".format(indx))
            eq.Residual = phi1_n[indx] - self.phiCC_n()
        for indx in range(N_p):
            eq = self.CreateEquation("phi1_p_{}".format(indx))
            eq.Residual = phi1_p[indx] - self.phiCC_p()

        # Set the solid average filling fraction to non-changing in the separator region.
        # The variable shouldn't actually be defined in the separator, but it's convenient
        # for plotting purposes that it be defined over the full domain, so we simply fix its value
        # in the separator.
        for indx in range(N_n, N_n + N_s):
            eq = self.CreateEquation("y_avg_s_{}".format(indx))
            eq.Residual = self.y_avg.dt(indx)

        # Define the total current.
        # There are multiple ways to do this. Here, we set the current to be the (negative of the)
        # integral of the reaction rate into all the particles in the negative electrode.
        # This is equivalent to setting it equal to
        # - the integral the reaction rate into all the particles in the positive electrode
        # - the current density in the electrolyte in the separator (which must be uniform)
        # - the current density in the solid bulk electrode at the current collector (if using
        #   finite conductivity in the electron-conducting phase)
        eq = self.CreateEquation("Total_Current")
        eq.Residual = self.current() + np.sum(self.F()*a[:N_n]*j_p[:N_n]*h_centers[:N_n])

        # Define the measured voltage
        # Thomas et al., below Eq 4
        eq = self.CreateEquation("Voltage")
        eq.Residual = self.V() - (self.phiCC_p() - self.phiCC_n())

        # For the simulation, we can either specify the voltage and let current be a calculated output,
        # or we can specify the current and let voltage be a calculated output.
        # These correspond to CV (constant voltage) and CC (constant current) operation respectively.
        # We ramp quickly from an equilibrium to the set point to facilitate the numerical calculation
        # of consistent initial conditions.
        if pinfo["profileType"] == "CC":
            # Total Current Constraint Equation
            eq = self.CreateEquation("Total_Current_Constraint")
            eq.Residual = self.current() - self.currset()*(1 - np.exp(-Time()/self.tau_ramp()))
        elif pinfo["profileType"] == "CV":
            # Keep applied potential constant
            eq = self.CreateEquation("applied_potential")
            eq.Residual = self.V() - self.Vset()*(1 - np.exp(-Time()/self.tau_ramp()))

class simTutorial(daeSimulation):
    def __init__(self):
        self.F = 96485.34 * A*s/mol
        # Define the model we're going to simulate
        # Constant current (CC) or constant voltage (CV) simulation
        profileType = "CC"
        # Time at which to stop the simulation (used for CV simulations)
        tend = 2*3200e0 * s
        # Fraction of battery capacity to (dis)charge (used for CC simulations)
        capfrac = 0.90
        # Applied current or voltage (used in CC and CV simulations respectively)
        self.currset = 3e+1 * A/m**2
        self.Vset = 3.5 * V
        # Lenghts of regions of battery, negative electrode, separator, positive electrode
        self.L_n = 243e-6 * m
        self.L_s = 50e-6 * m
        self.L_p = 200e-6 * m
        # Porosies in each region
        self.poros_n = 0.3
        self.poros_s = 1.0
        self.poros_p = 0.3
        self.L_tot = self.L_n + self.L_s + self.L_p
        # Number of grid points in each region
        self.N_n = 20
        self.N_s = 20
        self.N_p = 20
        # Number of radial grid points for active particles in each electrode
        self.NR_n = 15
        self.NR_p = 20
        # Radius of active particles in each electrode
        self.R_n = 18e-6 * m
        self.R_p = 1e-6 * m
        # Maximum concentration of Li in the active materials
        self.csmax_n = 13.2e3 * mol/m**3
        self.csmax_p = 23.72e3 * mol/m**3
        # Initial filling fraction of each electrode.
        # For discharge, negative starts full and positive empty, opposite for charge
        self.ff0_n = 0.99
        self.ff0_p = 0.21
        # Capacity of each electrode
        capacity_n = self.csmax_n*(1 - self.poros_n)*self.L_n*self.F
        capacity_p = self.csmax_p*(1 - self.poros_p)*self.L_p*self.F
        # Capacity of battery
        capacity = min(capacity_n, capacity_p)
        if profileType == "CC":
            tend = capfrac * capacity / self.currset
        self.process_info = {"profileType": profileType, "tend": tend}

        self.process_info["N_n"] = self.N_n
        self.process_info["N_s"] = self.N_s
        self.process_info["N_p"] = self.N_p
        self.m = ModCell("tutorial_che_6", process_info=self.process_info)
        self.m.Description = __doc__

    def SetUpParametersAndDomains(self):
        h_n = self.L_n / self.N_n
        h_s = self.L_s / self.N_s
        h_p = self.L_p / self.N_p
        xvec_centers_n = [h_n*(0.5 + indx) for indx in range(self.N_n)]
        xvec_centers_s = [self.L_n + h_s*(0.5 + indx) for indx in range(self.N_s)]
        xvec_centers_p = [(self.L_n + self.L_s) + h_p*(0.5 + indx) for indx in range(self.N_p)]
        xvec_centers = xvec_centers_n + xvec_centers_s + xvec_centers_p
        xvec_faces = [0 * m] + [h_n*(1 + indx) for indx in range(self.N_n)]
        xvec_faces += [self.L_n + h_s*(1 + indx) for indx in range(self.N_s)]
        xvec_faces += [(self.L_n + self.L_s) + h_p*(1 + indx) for indx in range(self.N_p)]
        # Domains in ModCell
        self.m.x_centers_n.CreateStructuredGrid(self.N_n - 1, 0, 1)
        self.m.x_centers_p.CreateStructuredGrid(self.N_p - 1, 0, 1)
        self.m.x_centers_full.CreateStructuredGrid(self.N_n + self.N_s + self.N_p - 1, 0, 1)
        self.m.x_faces_full.CreateStructuredGrid(self.N_n + self.N_s + self.N_p, 0, 1)
        self.m.x_centers_n.Points = [x.value for x in xvec_centers_n]
        self.m.x_centers_p.Points = [x.value for x in xvec_centers_p]
        self.m.x_centers_full.Points = [x.value for x in xvec_centers]
        self.m.x_faces_full.Points = [x.value for x in xvec_faces]
        # Domains in each particle
        for indx_n in range(self.m.x_centers_n.NumberOfPoints):
            self.m.particles_n[indx_n].r.CreateStructuredGrid(self.NR_n - 1, 0, self.R_n.value)
        for indx_p in range(self.m.x_centers_p.NumberOfPoints):
            self.m.particles_p[indx_p].r.CreateStructuredGrid(self.NR_p - 1, 0, self.R_p.value)
        # Parameters in ModCell
        self.m.R.SetValue(8.31447 * J/(mol*K))
        self.m.T.SetValue(298 * K)
        self.m.c_ref.SetValue(1000 * mol/m**3)
        self.m.tau_ramp.SetValue(1e-3 * self.process_info["tend"])
        # Parameters in each particle
        for indx_n in range(self.m.x_centers_n.NumberOfPoints):
            p = self.m.particles_n[indx_n]
            p.j_0.SetValue(1e-4 * mol/(m**2 * s))
        for indx_p in range(self.m.x_centers_p.NumberOfPoints):
            p = self.m.particles_p[indx_p]
            p.j_0.SetValue(1e-4 * mol/(m**2 * s))

    def SetUpVariables(self):
        cs0_n = self.ff0_n * self.csmax_n
        cs0_p = self.ff0_p * self.csmax_p
        # ModCell
        for indx in range(self.m.x_centers_full.NumberOfPoints):
            self.m.c.SetInitialCondition(indx, 1e3 * mol/m**3)
            if indx < self.N_n:
                self.m.y_avg.SetInitialCondition(indx, self.ff0_n)
            elif indx < self.N_n + self.N_s:
                self.m.y_avg.SetInitialCondition(indx, 0.)
            elif indx < self.N_n + self.N_s + self.N_p:
                self.m.y_avg.SetInitialCondition(indx, self.ff0_p)
        # particles
        for indx_n in range(self.m.x_centers_n.NumberOfPoints):
            p = self.m.particles_n[indx_n]
            for indx_r in range(1, p.r.NumberOfPoints-1):
                p.c.SetInitialCondition(indx_r, cs0_n)
        for indx_p in range(self.m.x_centers_p.NumberOfPoints):
            p = self.m.particles_p[indx_p]
            for indx_r in range(1, p.r.NumberOfPoints-1):
                p.c.SetInitialCondition(indx_r, cs0_p)

def run(**kwargs):
    # An exception will be raised at some point when the log function is called for a negative value
    simulation = simTutorial()
    timeHorizon = min(3770.0, simulation.process_info["tend"].value)
    return daeActivity.simulate(simulation, reportingInterval = timeHorizon / 100,
                                            timeHorizon       = timeHorizon,

if __name__ == "__main__":
    guiRun = False if (len(sys.argv) > 1 and sys.argv[1] == 'console') else True
    run(guiRun = guiRun)