尝试求解大型 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]

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

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)

抱歉，代码太长了，但我真的无法用其他方式显示它。无论如何，这是错误 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

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

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)

通过这些更改，代码运行s。它没有做任何有趣的事情，但这不足为奇，因为在这个阶段方式太复杂了。 18个方程只会混淆事物。我强烈建议你用 <= 2 个方程解决这个问题。

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

将来，请提供实际运行的示例（无论如何，直到错误点）。

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

Errors when trying to solve large system of PDEs

integration

pde

python-3.x

fipy