Search code examples
pythonfipy

Solving PDE with sources in Python

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:

enter image description here

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!

Solution

  • Note, as mentioned by wd15 in a comment, there is more in depth discussion and followup on the mailing list.

    First, both the initial conditions and the sources can all be made using the where logic.

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

    This avoids the explicit construction of the arrays.

    Second, there is a key difference between initial conditions and sources. In the original formulation, where you use the SetValue method on the field variable phi, you are providing an initial condition for the transient solution (or an initial guess for a direct steady state solve) instead of an actual source. So the correct approach is your second one in which you actually add the source/sink terms to the equation directly:

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

    Also, to attempt a direct steady solve, you would instead write the equation without the TransientTerm,

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

    then solve using eq.solve(), which will solve the combined DiffusionTerm and the sources/sinks.

    However, the direct-to-steady approach merits a few words of caution. First, there isn't really a good numerical way of directly computing steady state solutions for general systems. Often, your best bet is actually to solve the transient equation by time stepping from some initial condition until you reach steady state, as that's actually probably the most robust algorithm for solving for steady state profiles. In some cases, you can even do that with one large time step, as in the section beginning "fully implicit solutions are not without their pitfalls" here. Second, are you sure that your system admits a steady state? You have no-flux boundary conditions (implied because you didn't specify any other boundary conditions), but you have internal sources/sinks. Unless those sources/sinks produce/consume the field variable at exactly the same rate, you will have net accumulation in the system. A simple example with R = constant, uniform, and non-zero:

    dc/dt = R

    with no flux boundary conditions is an equation that doesn't admit any steady state. The addition of a diffusion term doesn't change that. If you were to add any dirichlet (specify value) boundary conditions, that would enable a steady solution, as the net production/consumption within the system could leave/enter via the system boundaries. A discussion on boundary conditions and how to apply them can be found here.

    Finally, another thing to note. The source/sink terms, as written are "zero-order," which will lead to, e.g. concentration going negative if the sink term is large enough and/or far enough away from the sources. If this occurs, the model clearly needs to be revised by, e.g., making the sink first-order,

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

    This would ensure that sinks "turn off" as phi goes to zero, but this could also by modified to couple to other field variables like local cell density.