Как связать PDE адвекционной диффузионной реакции с FiPy

Question

Я пытался решить 1D-связанные PDE для задачи адвекции-диффузии-реакции с помощью функции Matlab Pdepe ( https://www.mathworks.com/help/matlab/ref/pdepe.html). Эта функция не работает должным образом в моем случае высокого адвекционного члена по сравнению с диффузионным членом. Поэтому я искал и нашел эту возможность использования библиотеки Python FiPy для решения моей системы PDE. Мои начальные условия u1=1 для 4*L/10

Мои связанные уравнения имеют следующую форму:

du1 / dt = d / dx (D1 * du1 / dx) + g * x * du1 / dx - mu1 * u1 / (K + u1) * u2

du2 / dt = d / dx (D2 * du2 / dx) + g * x * du2 / dx + mu2 * u1 / (K + u1) * u2

Я попытался написать это, объединив примеры FiPy (examples.convection.exponential1DSource.mesh1D, examples.levelSet.advection.mesh1D, examples.cahnHilliard.mesh2DCoupled).

Следующие строки - не рабочий пример, а моя первая попытка написать код. Это мое первое использование FiPy (из тестов и примеров документации), так что это может показаться упущением для обычных пользователей.

from fipy import *

g = 0.66
L = 10.
nx = 1000
mu1 = 1.
mu2 = 1.
K = 1.
D1 = 1.
D2 = 1.

mesh = Grid1D(dx=L / 1000, nx=nx)

x = mesh.cellCenters[0]
convCoeff = g*(x-L/2)

u10 = 4*L/10 < x < 6*L/10
u20 = 1.

u1 = CellVariable(name="u1", mesh=mesh, value=u10)
u2 = CellVariable(name="u2", mesh=mesh, value=u20)

## Neumann boundary conditions
u1.faceGrad.constrain(0., where=mesh.facesLeft)
u1.faceGrad.constrain(0., where=mesh.facesRight)
u2.faceGrad.constrain(0., where=mesh.facesLeft)
u2.faceGrad.constrain(0., where=mesh.facesRight)

sourceCoeff1 = -1*mu1*u1/(K+u1)*u2
sourceCoeff2 = 1*mu2*u1/(K+u1)*u2

eq11 = (TransientTerm(var=u1) == DiffusionTerm(coeff=D1, var=u1) + ConvectionTerm(coeff=convCoeff))
eq21 = (TransientTerm(var=u2) == DiffusionTerm(coeff=D2, var=u2) + ConvectionTerm(coeff=convCoeff))

eq12 = ImplicitSourceTerm(coeff=sourceCoeff1, var=u1)
eq22 = ImplicitSourceTerm(coeff=sourceCoeff2, var=u2)

eq1 = eq11 & eq12
eq2 = eq21 & eq22

eqn = eq1 & eq2
vi = Viewer((u1, u2))

for t in range(100):
    u1.updateOld()
    u2.updateOld()
    eqn.solve(dt=1.e-3)
    vi.plot()

Спасибо за любое предложение или исправление. Если вам случится знать хорошее руководство по этой конкретной проблеме, я был бы рад прочитать его, так как я не нашел ничего лучше, чем примеры в документации FiPy.

Answer 1

Несколько вопросов:

• Цепные сравнения python не работают в numpy и, следовательно, не работают в FiPy. Итак, напишите

  u10 = (4*L/10 < x) & (x < 6*L/10)
  

  Кроме того, это делает u10 поле Booleans, которое смущает FiPy, так что пишите

  u10 = ((4*L/10 < x) & (x < 6*L/10)) * 1.
  

  или еще лучше написать

  u1 = CellVariable(name="u1", mesh=mesh, value=0., hasOld=True)
  u2 = CellVariable(name="u2", mesh=mesh, value=1., hasOld=True)
  u1.setValue(1., where=(4*L/10 < x) & (x < 6*L/10))
  

• ConvectionTerm принимает векторный коэффициент. Один из способов получить это

  convCoeff = g*(x-L/2) * [[1.]]
  

  которая представляет переменную 1D ранга 1
• Если вы объявите, какие Variable Term относится, вы должны сделать это для всех Term s, так напишите, например,

  ConvectionTerm(coeff=convCoeff, var=u1)
  

• ConvectionTerm(coeff=g*x, var=u1) не представляет g * x * du1/dx. Он представляет собой d(g * x * u1)/dx. Итак, я верю, что вы захотите

  ConvectionTerm(coeff=convCoeff, var=u1) - ImplicitSourceTerm(coeff=g, var=u1)
  

• ImplicitSourceTerm(coeff=sourceCoeff1, var=u1 не представляет -1*mu1*u1/(K+u1)*u2 скорее это представляет -1*mu1*u1/(K+u1)*u2*u1, Итак, для лучшей связи между уравнениями, напишите

  sourceCoeff1 = -mu1*u1/(K+u1)
  sourceCoeff2 = mu2*u2/(K+u1)
  
  ... ImplicitSourceTerm(coeff=sourceCoeff1, var=u2) ...
  ... ImplicitSourceTerm(coeff=sourceCoeff2, var=u1) ...
  

• Как указано @wd15 в комментариях, вы объявляете четыре уравнения для двух неизвестных. & не означает "сложить два уравнения вместе" (что можно сделать с помощью +), скорее это означает "решить эти два уравнения одновременно". Итак, напишите

  sourceCoeff1 = mu1*u1/(K+u1)
  sourceCoeff2 = mu2*u2/(K+u1)
  
  eq1 = (TransientTerm(var=u1) 
         == DiffusionTerm(coeff=D1, var=u1) 
         + ConvectionTerm(coeff=convCoeff, var=u1) 
         - ImplicitSourceTerm(coeff=g, var=u1) 
         - ImplicitSourceTerm(coeff=sourceCoeff1, var=u2))
  eq2 = (TransientTerm(var=u2) 
         == DiffusionTerm(coeff=D2, var=u2) 
         + ConvectionTerm(coeff=convCoeff, var=u2) 
         - ImplicitSourceTerm(coeff=g, var=u2) 
         + ImplicitSourceTerm(coeff=sourceCoeff2, var=u1))
  
  eqn = eq1 & eq2
  

• CellVariable должен быть объявлен с hasOld=True для того, чтобы позвонить updateOld(), так

  u1 = CellVariable(name="u1", mesh=mesh, value=u10, hasOld=True)
  u2 = CellVariable(name="u2", mesh=mesh, value=u20, hasOld=True)
  

Полный код, который, кажется, работает