Ошибки при попытке решить большую систему PDE - PullRequest
Новая
       65

Ошибки при попытке решить большую систему PDE

1 голос
/ 27 апреля 2019

У меня есть система из 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)

Извините за длинный код, но я не могу показать его по-другому.В любом случае это ошибка, которую он возвращает:

Файл "C: \ Users \ tjcze \ Anaconda3 \ lib \ site-packages \ fipy \ matrices \ scipyMatrix.py", строка 218, в addAt assert (len (id1) == len (id2) == len (vector))

AssertionError

Что я должен сделать, чтобы эта система работала правильно?

1 Ответ

0 голосов
/ 27 апреля 2019

Ошибка из-за того, что вы пытаетесь сделать с u_x.u_x должен иметь ранг-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 уравнений только запутывают вещи.Я настоятельно рекомендую устранить эту проблему с помощью уравнений <= 2. </p>

На этом этапе ваши уравнения вообще не связаны (eq0 только неявно зависит от EDC, eq1 только на EC и т. Д.).Это не так, но не очень полезно.Конечно, нет смысла в синтаксисе eq = eq0 & eq1 & ....EDC - единственная переменная, которая может эволюционировать, и она постоянна, поэтому она тоже не будет развиваться.

В будущем, пожалуйста, предоставьте примеры, которые действительно выполняются (в любом случае, вплоть до ошибки).

...