PythonソースでPDEを解決する

Question

私はBiologyに触発された問題に対処するために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 

ソースおよびシンクが提供する「無限」ソースとシンクこれを変更することはありませんことを考えます。

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 

を使用して、定常状態を解くしようとすると、私が取得します。しかし、私は再び使用「の手順で」それを解決した場合：私は期待しているはずだ何

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

が、私はと類似の素敵な画像を取得します：

私が指定する方法上の任意のアドバイス定常状態の解を得ることができるような方法で異なる放射/消費率を持つ、異なる空間的位置にあるソース/シンクとの初期設定？

ありがとうございます！

Answer 1

コメントでwd15で言及されているように、メーリングリストについての詳細な議論とフォローアップがあります。

最初に、初期条件とソースの両方をロジックwhereで作成することができます。

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

これにより、配列が明示的に構築されるのを回避します。

第2に、初期条件とソースには大きな違いがあります。フィールド変数phiにSetValueメソッドを使用する元の形式では、実際のソースではなく、過渡解の初期条件（または直接定常状態解の初期推測）を提供しています。だから、正しいアプローチは、あなたが実際に直接方程式にソース/シンク用語を追加した2番目のいずれかになります。

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

また、解決直接安定を図るために、あなたの代わりにTransientTermせずに式を記述し、

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

eq.solve()を使用して解決すると、DiffusionTermとソース/シンクの組み合わせが解決します。

しかし、直接的に安定したアプローチには注意が必要です。まず、一般的なシステムの定常解を直接計算するための数値的な方法は実際にはありません。多くの場合、実際には定常状態に達するまで、ある初期状態から時間をずらして過渡方程式を解くことが最善の方法です。場合によっては、「暗黙の完全なソリューションに落とし穴がないわけではありません」というセクションのように、1つの大きな時間ステップでこれを実行することもできます。here。第二に、あなたのシステムが安定した状態になっていることを確信していますか？あなたはフラックス境界条件を持っていません（あなたが他の境界条件を指定していないために暗示されます）が、内部にソース/シンクがあります。これらのソース/シンクがまったく同じレートでフィールド変数を生成/消費しない限り、システムには正味の蓄積があります。 R =定数、均一な、非ゼロの単純な例：

dc/dt = R

ないフラックス境界条件は、任意の定常状態を認めない式です。拡散項の追加はそれを変えない。任意のディリクレ（値を指定する）境界条件を追加すると、システム内の正味の生産/消費がシステム境界を離れて出入りすることがあるため、安定した解決策が可能になります。境界条件とその適用方法については、hereを参照してください。

最後に、別のことに注意してください。書かれているように、ソース/シンク用語は「ゼロ次」であり、これは例えば次のようになる。シンクの項が十分に大きい場合、および/またはソースから十分に離れている場合、濃度は負になる。この問題が発生した場合、モデルは明らかに

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

これはPHIがゼロになるよう、「オフ」沈むという保証するであろう、シンク一次を作り、例えば、によって改訂される必要があり、これもによって変更できローカルセル密度のような他のフィールド変数に結合することができます。