Решение PDE с исходниками в Python

Я использую FiPy для решения проблемы, вдохновленной биологией.

По сути, я хочу представить 2D-плоскость, где в разных точках у меня есть источники и стоки. Источники выделяют субстрат с фиксированной скоростью (у разных источников может быть разная скорость), а поглотители потребляют субстрат с фиксированной скоростью (у разных поглотителей может быть разная скорость). Мой код:

import numpy.matlib
from fipy import CellVariable, Grid2D, Viewer, TransientTerm, DiffusionTerm, ImplicitSourceTerm, ExplicitDiffusionTerm
from fipy.tools import numerix
from time import *

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

arr_grid = numerix.array((
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,),'d')

arr_source = numerix.array((
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0.5,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,),'d')

arr_sink = numerix.array((
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0.5,),'d')

source = CellVariable(mesh=mesh, value = arr_source)
sink = CellVariable(mesh=mesh, value = arr_sink)
phi = CellVariable(name = "solution variable", mesh = mesh, value = arr_grid)
X,Y = mesh.cellCenters
phi.setValue(3.0, where=(X < 2.0) & (X > 1.0))
phi.setValue(-1.0, where=(X < 6.0) & (X > 5.0))
D = 1.
eq = TransientTerm() == DiffusionTerm(coeff=D)



viewer = Viewer(vars=phi, datamin=0., datamax=1.)

steadyState = False

if(steadyState):
    print("SteadyState")
    DiffusionTerm().solve(var=phi)
    viewer.plot()
    sleep(20)
else:
    print("ByStep")
    timeStepDuration = 10 * 0.9 * dx**2 / (2 * D)
    steps = 500
    for step in range(steps):
        print(step)
        eq.solve(var=phi,
        dt=timeStepDuration)
        if __name__ == '__main__':
            viewer.plot()

Это работает хорошо, но FiPy рассматривает источники как «невозобновляемые», и в конечном итоге я получаю однородную концентрацию по всему пространству, как и следовало ожидать. Альтернативой было удаление:

X,Y = mesh.cellCenters
phi.setValue(3.0, where=(X < 2.0) & (X > 1.0))
phi.setValue(-1.0, where=(X < 6.0) & (X > 5.0))

И измените уравнение на:

eq = TransientTerm() == DiffusionTerm(coeff=D) + source - sink

Учитывая, что источник и приемник никогда не меняются, предлагаются «бесконечные» источники и приемники. Однако, когда я пытаюсь решить для установившегося состояния, используя

eq = TransientTerm() == DiffusionTerm(coeff=D) + source - sink

Я получил:

C:\Python27\python.exe C:/Users/dario_000/PycharmProjects/mesa-test/mesher.py
SteadyState
C:\Python27\lib\site-packages\fipy-3.1.dev134+g64f7866-py2.7.egg\fipy\solvers\scipy\linearLUSolver.py:71: RuntimeWarning: invalid value encountered in double_scalars
  if (numerix.sqrt(numerix.sum(errorVector**2)) / error0)  <= self.tolerance:

И уравнение не решается. Однако, если я решу это «по шагам», используя снова:

eq = TransientTerm() == DiffusionTerm(coeff=D) + source - sink

Я получаю красивую картинку, похожую на то, что я ожидал:

введите здесь описание изображения

Любые советы о том, как я могу указать первоначальную настройку с источниками/стоками в разных пространственных положениях, каждый из которых имеет разную скорость выбросов/потребления таким образом, чтобы я мог получить стационарное решение?

Спасибо!


python fipy
person MrD    schedule 03.04.2016    source источник
comment
Этот вопрос дублируется в списке рассылки FiPy, где есть больше обсуждений, см. ветку . gmane.org/gmane.comp.python.fipy/4034.   -  person wd15    schedule 04.04.2016

Ответы (1)


arrow_upward
4
arrow_downward

Обратите внимание, как упоминалось в комментарии wd15, в списке рассылки есть более подробное обсуждение и последующие действия.

Во-первых, и начальные условия, и источники могут быть созданы с использованием логики where.

source = CellVariable(mesh=mesh, value = arr_source, where=(2 < X) & (X < 3))

Это позволяет избежать явного построения массивов.

Во-вторых, существует ключевое различие между начальными условиями и источниками. В исходной формулировке, где вы используете метод SetValue для переменной поля phi, вы предоставляете начальное условие для переходного решения (или начальное предположение для прямого решения в установившемся режиме) вместо фактического источника. Таким образом, правильный подход - ваш второй, в котором вы фактически добавляете условия источника/приемника к уравнению напрямую:

eq = TransientTerm() == DiffusionTerm(coeff=D) + source - sink

Кроме того, чтобы попытаться решить прямое устойчивое решение, вместо этого вы должны написать уравнение без TransientTerm,

eq = 0 == DiffusionTerm(coeff=D) + source - sink

затем решите, используя eq.solve(), который решит комбинированный DiffusionTerm и источники/приемники.

Тем не менее, подход «прямой к устойчивому» заслуживает нескольких предостережений. Во-первых, на самом деле не существует хорошего численного способа прямого вычисления стационарных решений для общих систем. Часто лучше всего решить уравнение переходного процесса путем шага по времени от некоторого начального состояния до достижения устойчивого состояния, поскольку на самом деле это, вероятно, самый надежный алгоритм для решения профилей установившегося состояния. В некоторых случаях вы даже можете сделать это с одним большим временным шагом, как в начале раздела «полностью неявные решения не лишены ловушек» здесь. Во-вторых, уверены ли вы, что ваша система допускает устойчивое состояние? У вас есть граничные условия отсутствия потока (подразумеваемые, потому что вы не указали никаких других граничных условий), но у вас есть внутренние источники/стоки. Если только эти источники/приемники не производят/не потребляют переменную поля с одинаковой скоростью, у вас будет чистое накопление в системе. Простой пример с R = постоянным, однородным и ненулевым:

dc/dt = R

без граничных условий потока - это уравнение, которое не допускает стационарного состояния. Добавление термина диффузии не меняет этого. Если бы вы добавили какие-либо граничные условия дирихле (укажите значение), это обеспечило бы устойчивое решение, поскольку чистое производство/потребление в системе могло бы выходить/входить через границы системы. Обсуждение граничных условий и способов их применения можно найти здесь< /а>.

Наконец, еще одно замечание. Условия источника/приемника, как написано, являются «нулевым порядком», что приведет, например, к концентрация становится отрицательной, если член стока достаточно велик и/или достаточно далеко от источников. Если это произойдет, то модель явно необходимо пересмотреть, например, сделав сток первого порядка,

eq = TransientTerm() == DiffusionTerm(coeff=D) + source - sink*phi

Это гарантировало бы, что стоки «выключаются», когда фи приближается к нулю, но это также может быть изменено для связи с другими переменными поля, такими как локальная плотность клеток.

person muon    schedule 05.04.2016