Решение PDE с источниками в Python

Question

Я использую 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 

Я получаю хорошую картину, сходную с тем, что я бы ожидать:

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

Спасибо!

Answer 1

Обратите внимание, что, как упоминалось в комментарии 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 и источники/поглотители.

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

dc/dt = R

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

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

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

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