Errors when trying to solve large system of PDEs

I have a system of 18 linear (i'm assuming) PDEs. Howe ever they are very awkwardly shaped PDEs. Every time I try to use FiPy to solve the system of coupled PDEs I get an error. I have fixed most of them on my own but I can't seem to get past this one.

First you will need a list of equations that will be added later in the program:

import time
from fipy import Variable, FaceVariable, CellVariable, Grid1D, Grid2D, ExplicitDiffusionTerm, TransientTerm, DiffusionTerm, Viewer,PowerLawConvectionTerm, ImplicitSourceTerm
from fipy.tools import numerix
import math
import numpy as np
import sympy as sp
Temperature = 500.0 #Temperature in Celsius

Temp_K = Temperature + 273.15 #Temperature in Kelvin

Ea = [342,7,42,45,34] #Activation energy in order from the list excel/word file
#Frequency factor ko (Initial k)
k_0 =[5.9 * (10**15), 1.3 * (10**13), 1.0 * (10**12), 5.0 * (10**11), 1.2 * (10**13)]
fk1 = [math.exp((-1.0 * x)/(R*Temp_K)) for x in Ea] #Determines k value at given temperature value (exp(-Ea/R*T))
final_k = [fk1*k_0 for fk1,k_0 in zip(fk1,k_0)] #Multiplys by the inital k value to determine the k value at the given temperature (ko * exp(-Ea/R*T))
     final_kcm = [x for x in final_k]
def rhs(eq):
       eq_list = [-1.0*final_kcm[0]*EDC - final_kcm[1]*EDC*R1 - final_kcm[2]*EDC*R2 - final_kcm[3]*EDC*R4 - final_kcm[4]*EDC*R5 - final_kcm[5]*EDC*R6,
       final_kcm[2]*EDC*R2 - final_kcm[8]*EC*R1 + final_kcm[13]*VCM*R2,
       final_kcm[1]*EDC*R1 + final_kcm[6]*R2*R1 + final_kcm[7]*R1*R3 + final_kcm[8]*R1*EC + final_kcm[9]*R1*C11 + final_kcm[10]*R1*C112 +  final_kcm[12]*R1*VCM + 2.0*final_kcm[20]*R1*C2H2,
       2.0*final_kcm[20]*R2*C2H2,
       final_kcm[15]*R5*VCM,

       return eq_list[eq]

Here is the rest of the program I use to generate the system of PDEs.

EDC = CellVariable(mesh=mesh, hasOld=True, value=10)
EC = CellVariable(mesh=mesh, hasOld=True, value=0)
HCl = CellVariable(mesh=mesh, hasOld=True, value=0)
Coke = CellVariable(mesh=mesh, hasOld=True, value=0)
CP = CellVariable(mesh=mesh, hasOld=True, value=0)

EDC.constrain(10, mesh.facesLeft)
EC.constrain(0., mesh.facesLeft)
HCl.constrain(0., mesh.facesLeft)
Coke.constrain(0., mesh.facesLeft)
CP.constrain(0., mesh.facesLeft)

nsp =18
u_x = [[ [0,]*nsp for n in range(nsp) ]]
for z in range(nsp):
    u_x[0][z][z] = 1.0

eq0 = TransientTerm(var = EDC) == PowerLawConvectionTerm(coeff = u_x, var = EDC) + ImplicitSourceTerm(rhs(0),var = EDC)
eq1 = TransientTerm(var = EC) == PowerLawConvectionTerm(coeff = u_x, var = EC) + ImplicitSourceTerm(rhs(1),var = (EC))
eq2 = TransientTerm(var = HCl) == PowerLawConvectionTerm(coeff = u_x, var = HCl) + ImplicitSourceTerm(rhs(2),var = (HCl))
eq3 = TransientTerm(var = Coke) == PowerLawConvectionTerm(coeff = u_x, var = Coke) + ImplicitSourceTerm(rhs(3),var = (Coke))
eq4 = TransientTerm(var = CP) == PowerLawConvectionTerm(coeff = u_x, var = CP) + ImplicitSourceTerm(rhs(4),var = (CP))


eqn = eq0 & eq1 & eq2 & eq3 & eq4 

if __name__ == '__main__':
     viewer = Viewer(vars = (EDC,EC,HCl,Coke,CP))
     viewer.plot()

for t in range(1):
    EDC.updateOld()
    EC.updateOld()
    HCl.updateOld()
    Coke.updateOld()
    CP.updateOld()

    eqn.solve(dt=1.e-3)

Sorry for the long code but I can't really show it any other way. Anyway this is the error it returns :

File "C:\Users\tjcze\Anaconda3\lib\site-packages\fipy\matrices\scipyMatrix.py", line 218, in addAt assert(len(id1) == len(id2) == len(vector))

AssertionError

What should I do to get this system to work correctly?

Solution

The error is because of whatever you're trying to do with u_x. u_x should be a rank-1 FaceVariable. Its shape should be #dims x #faces (or #dims x 1); it should not be (#eqns x #eqns).

Setting u_x = [1.] gets rid of the AssertionError.

You will then get a series of warnings:

UserWarning: sweep() or solve() are likely to produce erroneous results when `var` does not contain floats.

Fix this by initializing all of your CellVariables with floats, e.g.,

EDC = CellVariable(mesh=mesh, hasOld=True, value=10.)

instead of

EDC = CellVariable(mesh=mesh, hasOld=True, value=10)

With those changes, the code runs. It doesn't do anything interesting, but that's hardly surprising, as it's way too complicated at this stage. 18 equations only obfuscates things. I strongly recommend you troubleshoot this problem with <= 2 equations.

At this point, your equations aren't coupled at all (eq0 only implicitly depends on EDC, eq1 only on EC, etc.). This isn't wrong, but not terribly useful. Certainly no point in the eq = eq0 & eq1 & ... syntax. EDC is the only variable with a chance of evolving, and it's constant, so it won't evolve, either.

In future, please provide examples that actually run (up to the point of the error, anyway).