Hi,
I originally asked this on
http://stackoverflow.com/questions/36385367/solving-pde-with-sources-in-python
but am copying it to this list in case somebody has any advice :) I've
already read through
https://www.mail-archive.com/[email protected]/msg02439.html but didn't solve
my problem.
I'm using FiPy to address a problem inspired by Biology.
Essentially, I want to represent a 2D plane where at different points I
have sources and sinks. Sources emit substrate at a fixed rate (different
sources can have different rates) and sinks consume substrate at a fixed
rate (different sinks can have different rates). My code:
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()
This works well, but FiPy treats the sources as "non-renewable" and
eventually I get an homogeneous concentration throughout the space as would
be expected. The alternative was to delete:
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))
And change the equation to:
eq = TransientTerm() == DiffusionTerm(coeff=D) + source - sink
Given that source and sink are never changed this offered "infinite"
sources and sinks. However, when I try solving for steady-state using
eq = TransientTerm() == DiffusionTerm(coeff=D) + source - sink
I get:
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:
And the equation isn't solved. However if I solve it "by steps" using again:
eq = TransientTerm() == DiffusionTerm(coeff=D) + source - sink
I get a nice picture similar to what I'd be expecting:
http://i.stack.imgur.com/OsJda.png
Any advice on how I can specify an initial setup with sources/sinks at
different spatial positions each with a different emission/consumption rate
in such a way that I may be able to obtain the steady-state solution?
Thanks!
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
