Блог

Я работаю с биологической моделью распределения микробной биомассы (b1) на 2D-сетке. Из биомассы образуется белок (р1). Биомасса диффундирует по сетке, в то время как белок этого не делает. Предполагается, что биомасса диффундирует только в том случае, если вырабатывается определенное количество белка (p > p_lim). Я пытаюсь реализовать это, используя фиктивную переменную ячейки z, умноженную на коэффициент диффузии, и устанавливая его от 0 до 1 только в ячейках, где p> p_lim.

Условие работает нормально, и когда в ячейке достигается критическое количество p, z устанавливается равным 1, и происходит диффузия. Однако диффузия по-прежнему не работает с той скоростью, которую я хотел бы, потому что для вычисления диффузии используется переменная face, а не значение самой ячейки. Грани z всегда являются средним значением ячейки с z = 1 и соседних ячеек с z = 0. I Я, однако, хотел бы, чтобы диффузия работала с первоначальной скоростью, даже если соседняя ячейка все еще находится на уровне p < p_lim.

Итак, мой вопрос таков: могу ли я каким-то образом получить доступ к faceVariable и изменить его? Например, установите для грани значение 1, если какая-либо соседняя ячейка достигла p1> p_lim? Я предполагаю, что это не совсем математическая вещь, но я не мог придумать другого способа смоделировать эту проблему.

Ниже я покажу очень уменьшенную форму моей модели. В любом случае, я вам очень благодарен за уделенное время!

 ##### produce mesh

nx= 5.
ny= nx
dx = 1.
dy = dx
L = nx*dx
mesh = Grid2D(nx=nx,ny=ny,dx=dx,dy=dy)


#parameters
h1 = 0.5 # production rate of p
Db = 10. # diffusion coeff of b
p_lim=0.1 


# cell variables
z = CellVariable(name="z",mesh=mesh,value=0.)

b1 = CellVariable(name="b1",mesh=mesh,hasOld=True,value=0.)

p1= CellVariable(name="p1",mesh=mesh,hasOld=True,value=0.)


# equations
eqb1 = (TransientTerm(var=b1)== DiffusionTerm(var=b1,coeff=Db*z.arithmeticFaceValue)-ImplicitSourceTerm(var=b1,coeff=h1))
eqp1 = (TransientTerm(var=p1)==ImplicitSourceTerm(var=b1,coeff=h1)) 

# set b1 to 10. in the center of the grid
b1.setValue(10.,where=((x>2.)amp;(x<3.)amp;(y>2.)amp;(y<3.)))
         
vi=Viewer(vars=(b1,p1),FIPY_VIEWER="matplotlib")


eq = eqb1 amp; eqp1

from builtins import range
for t in range(10):
    b1.updateOld()
    p1.updateOld()
    z.setValue(z   0.1,where=((p1>=p_lim) amp; (z < 1.)))
    
    eq.solve(dt=0.1)
     
    vi.plot()
 

В дополнение к .arithmeticFaceValue , FiPy предоставляет другие интерполяторы между значениями ячейки и номинала, такие как .harmonicFaceValue и .minmodFaceValue .

Эти свойства реализованы с использованием подклассов _CellToFaceVariable , в частности _ArithmeticCellToFaceVariable , _HarmonicCellToFaceVariable , и _MinmodCellToFaceVariable .

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

Вам необходимо переопределить _calc_() метод, чтобы предоставить свой пользовательский расчет. Этот метод принимает три аргумента:

• alpha : массив отношения (0-1) расстояния от грани до ячейки с одной стороны относительно расстояния от расстояния от ячейки с другой стороны до ячейки с первой стороны
• id1 : массив индексов ячеек на одной стороне грани
• id2 : массив индексов ячеек на другой стороне грани

Примечание: вы можете игнорировать любое предложение if inline.doInline: и посмотреть на _calc_() метод, определенный в else: предложении.