尝试求解大型 PDE 系统时出错

Question

我有一个由 18 个线性（我假设）PDE 组成的系统。但是，它们是形状非常笨拙的 PDE。每次我尝试使用 FiPy 解决耦合 PDE 的系统时，都会出现错误。我自己修复了其中的大部分，但我似乎无法克服这个问题。

首先，您需要一个方程式列表，稍后将在程序中添加：

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]

这是我用来生成偏微分方程系统的程序的其余部分。

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)

抱歉，代码太长了，但我真的无法以其他方式展示它。无论如何，这是它返回的错误：

文件 "C:\Users\tjcze\Anaconda3\lib\site-packages\fipy\matrices\scipyMatrix.py", 第 218 行，在 addAt 断言(len(id1) == len(id2) == len(vector))

断言错误

我应该怎么做才能让这个系统正常工作？

Answer 1

该错误是因为您尝试对u_x 执行任何操作。 u_x 应该是 rank-1 FaceVariable。 它的形状应该是#dims x #faces（或#dims x 1）；它不应该是（#eqns x #eqns）。

设置u_x = [1.] 去掉了AssertionError。

然后您将收到一系列警告：

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

通过使用浮点数初始化所有 CellVariables 来解决此问题，例如，

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

而不是

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

通过这些更改，代码运行。它没有做任何有趣的事情，但这并不奇怪，因为在这个阶段它方式太复杂了。 18 个方程只会混淆事物。我强烈建议您使用

此时，您的方程根本没有耦合（eq0 仅隐式依赖于 EDC，eq1 仅依赖于 EC 等）。这没有错，但不是非常有用。当然eq = eq0 & eq1 & ... 语法毫无意义。 EDC 是唯一有机会进化的变量，而且它是恒定的，所以它也不会进化。

以后，请提供实际运行的示例（无论如何，直到出现错误为止）。