FiPy system of coupled non-linear PDEs - Error in convection coefficients

Question

Edit: Problem is being discussed on Github(https://github.com/usnistgov/fipy/issues/883)

I am currently learning how to model in FiPy. I am trying to reproduce a model from a paper on oxygen evolution in porous electrodes. I keep getting the error that the convection terms require a vector as a coefficient. Also, I am not sure if such a complex coupled model can be sovled with FiPy.

The system of differential equations describes the formation and transport of oxygen in a porous electrode in 1-D.

The code is below, here are links to the equations (reputation not high enough to post as pictures directly). I want to solve for the electrolyte concentration c, the oxygen concentration cO2, the current i2 and the electrolyte potential phi2.

Equations

Here is how I tried to implement it in FiPy (: Edit: cleaned up based on a nice example I found here:

#  %%
from fipy import *
from fipy import CellVariable, Variable, Grid1D, TransientTerm, DiffusionTerm, ImplicitSourceTerm, LinearLUSolver, Viewer
import math
from fipy.tools import numerix
from fipy.solvers.scipy import LinearPCGSolver, LinearLUSolver
import numpy as np

#%%
# Define constants 

DOH = 1e-3  # [cm2/s] # test
DO2 = 1e-3 # [cm2/s]
t0 = 0.78
eps = 0.4
F = 96485.3321
R = 8.314
T = 298
gamma = 1.5
kappa = 12 
cOH_0 = 7.1e-3
a = 3864 # [cm2/cm3]
io2ref = 1e-11
alpha_an = 1.5
alpha_cat = 0.5
cO2 = 1e-7
cO2ref = 1e-7
phi1 = 0 # arbitrary value 
Uref = 0.3027-0.0983
iapp = 0.200 # A/cm2

# %%
# Mesh
L = 1e-3
nx = 1000
msh = Grid1D(dx=L / nx, nx=nx)
x = msh.cellCenters[0]

# initial conditions

cOH = CellVariable(name='OH concentration', mesh=msh, value= 7.0e-3, hasOld= True)
cO2 = CellVariable(name='O2 concentration', mesh=msh, value= 0., hasOld=True)
i2 = CellVariable(name='pore current', mesh=msh, value= 0.,hasOld=True)
phi2 = CellVariable(name='electrolyte potential', mesh=msh, value= 0.,hasOld=True)

# Boundary Conditions
i2.constrain(0, where=msh.facesLeft)
i2.constrain(iapp, where=msh.facesRight)
cOH.constrain(cOH_0, where=msh.facesRight)
cO2.constrain(0, where=msh.facesRight)
phi2.faceGrad.constrain(-iapp/(kappa*eps**gamma), msh.facesRight)

# Source

j2 = a*io2ref*((cOH/cOH_0)**2*numerix.exp(alpha_an*F*(phi1-phi2-Uref )/(R*T))-(cO2/cO2ref)*numerix.exp(-alpha_cat*F*(phi1-phi2-Uref )))
coeff_Ohm = -2*R*T/(cOH*F)*(1-t0+cOH/(2*cOH_0))
DO2_eff = eps**gamma * DO2
DOH_eff = eps**gamma*DOH
kappa_eff = kappa*eps**gamma

# Equations

eq1 = (TransientTerm(var=cOH, coeff=eps) == DiffusionTerm(var=cOH, coeff = DOH_eff) + (t0-1)/F*cOH.grad)

eq2 =(TransientTerm(var=cO2, coeff = eps) == DiffusionTerm(var=cO2, coeff=DO2_eff) + 1/(4*F)*j2)

eq3 =  i2/(kappa_eff) == ConvectionTerm(var = phi2, coeff = [[-1.]]) + ConvectionTerm(var = cOH, coeff = [[coeff_Ohm]]) - ImplicitSourceTerm(var = cOH, coeff=coeff_Ohm.grad)

eq4 = (i2.grad == j2)

#%%
vi = Viewer((cOH))

dt = 1e-3

eqsum = eq1 & eq2  & eq3 & eq4

solver = LinearLUSolver()
for t in range(100):
    cOH.updateOld()
    cO2.updateOld()
    phi2.updateOld()
    i2.updateOld()
    
    for sweep in range(5):
        
    eqsum.sweep(dt=dt,solver=solver)
        
        print('res1', res1)
        print('res2', res2)
        print('res3', res3)
        print('res3', res4)
    vi.plot()
# %%

I am not sure when I need to express the simple gradient of my solution variabel as a convection term. For example in equation 3. Should I write is as coeff * cOH.grad or ConvectionTerm(var=cOH, coeff = coeff)? Even though I write the coefficients as vectors using double square brackets, I get the error message that it needs to be a vector.

If I write the gradients in x simply as cOH.grad, then I get the following error:

TypeError: The coefficient must be rank 0 for a rank 0 solution variable.

I learned that the coefficient in the term ConvectioTerm needs to be in the derivative, so d(coeff*coH)/dx. In equation 3 the factor in front of dcOH/dx, let's call it f, also depends on my solution variable. I saw in this post that this can be written in the form of a convection term and a source term:

ConvectionTerm(var=cOH, coeff = f) - ImplicitSourceTerm(var=cOH, coeff = f.grad)

But when i want to couple the equations I get the following error:

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

Thanks a lot for your help and for taking the time to read this and going through it.