Steady-Source solution of PDE with fixed sources and sinks

[Dario Panada]

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/fi...@list.nist.gov/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!

[Raymond Smith]

It seems there's a bit of confusion here about initial condition vs. source
term.
Some comments in line.
Hope it helps.

Ray

On Sun, Apr 3, 2016 at 8:38 AM, Dario Panada <dario.panada at gmail.com> wrote:

> 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/fi...@list.nist.gov/msg02439.html but didn't solve
> my problem.
>

The list archives can be a good source, and I'd probably even more highly
recommend the examples on the fipy webpage in terms of learning per time.

As arr_grid is simply zeros, you could write
phi = CellVariable(name="solution variable", mesh=mesh, value=0.)
and avoid arr_grid entirely.

> 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))
>

You might consider setting up your source/sink CellVariables with this same
sort of "where" logic so that you avoid writing full arrays out explicitly.
What you've done with this line is establish an initial condition for the
field variable (or an initial guess if you're trying to do a steady-state
solve directly). Also, because this is an initial condition -- assuming phi
is representing, e.g., concentration -- then, it's strange that phi is set
to negative. I think this gets back to the issue that this is an initial
condition and not a source/sink term.

> D = 1.
> eq = TransientTerm() == DiffusionTerm(coeff=D)
>

There is no source/sink in this equation, so the steady profile will be a
perfectly uniform system with the average system concentration defined by
the initial conditions. Actually, I'm not sure how FiPy treats the
steady-state initial guess for Laplace's equation with no flux boundary
conditions like yours here. The governing equation + BC's without an
initial condition admits any uniform profile as a solution. Anyway, the
approach you took below does have the source/sink included in the actual
equations, so that seems like the direction to pursue.

>
>
>
> 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:
>

Here, I'm not sure what you mean by "non-renewable." As written above,
there isn't any source/sink at all.

>
> 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 don't see any difference here. Did you forget to edit this line?
To try to solve directly for steady state, you could write
eq = DiffusionTerm(coeff=D) + source - sink
then run
eq.solve()
noting that you call the eq solve method, not one from only the
DiffusionTerm as written above. I'm still unsure about the treatment of the
initial phi values in this sense as mentioned above.
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. Second, are you sure that
your system admits a steady state? You have no-flux 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. Simple example with R = constant and
uniform, non-zero value:
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.

>
> 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
> fi...@list.nist.gov
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://groups.google.com/a/list.nist.gov/forum/#!forum/fipy ]
>
>

[Dario Panada]

Hi Ray,

Many thanks for your e-mail and for your suggestions, I have made some
changes to your code as you suggested:

Firstly, I have re-written my CellVariables so I don't have to write the
array in full each time:

source = CellVariable(mesh=mesh, value=0.)
sink = CellVariable(mesh=mesh, value=0.)

phi = CellVariable(name="solution variable", mesh=mesh, value=0.)

X,Y = mesh.cellCenters
source.setValue(100.0, where=(X < 2.0) & (X > 1.0) & (Y > 1.0) & (Y < 2.0))

By "renewable" I meant that at each iteration the value at the source-point
remains the same (Ie: Acting as an infinite source emitting substrate at a
fixed rate, I hope this makes sense...?), I think I achieved this by using
the source CellVariable object which doesn't change.

I then added a no-flux constraint along the edges:

phi.constrain(0., where=mesh.exteriorFaces)

I am using the equation:

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

I then used the iterative approach as you suggested, and for this setup
obtain the following after about 60 iterations:

[image: Inline image 1]

Which is a bit surprising. As I only have a source I was expecting the
whole environment to "flood" with substrate but maybe this is just due to
me having misunderstood how this PDE is working?

If I then add a sink as

sink.setValue(50.0, where=(X > 8.0) & (X < 9.0) & (Y > 8.0) & (Y < 9.0))

I obtain the following after about 30 iterations:

[image: Inline image 2]

My full code: http://pastebin.com/8HUPh7xx

Does this look reasonable or do you think I should still fix something...?

Kind Regards,
Dario

[Raymond Smith]

Hi, Dario.

A few notes in line.

On Sun, Apr 3, 2016 at 2:17 PM, Dario Panada <dario.panada at gmail.com> wrote:

> Hi Ray,
>
> Many thanks for your e-mail and for your suggestions, I have made some
> changes to your code as you suggested:
>
> Firstly, I have re-written my CellVariables so I don't have to write the
> array in full each time:
>
> source = CellVariable(mesh=mesh, value=0.)
> sink = CellVariable(mesh=mesh, value=0.)
> phi = CellVariable(name="solution variable", mesh=mesh, value=0.)
> X,Y = mesh.cellCenters
> source.setValue(100.0, where=(X < 2.0) & (X > 1.0) & (Y > 1.0) & (Y < 2.0))
>
> By "renewable" I meant that at each iteration the value at the
> source-point remains the same (Ie: Acting as an infinite source emitting
> substrate at a fixed rate, I hope this makes sense...?), I think I achieved
> this by using the source CellVariable object which doesn't change.
>

This is contradictory. With a constant source term, the _value_ of phi
should _not_ remain constant in time at the source point. See some comments
below to help see this in your own simulations. An infinite source pumps
stuff into the system locally at a fixed rate, which is very different from
holding the value of phi fixed at a point. If you're actually trying to
model a source, then you're on the right track as-is. If you're looking to
constrain phi to have a fixed value internally, then the current
implementation will not do that. To do that, start by reading over the
"applying internal boundary conditions" section here:
http://www.ctcms.nist.gov/fipy/documentation/USAGE.html

>
> I then added a no-flux constraint along the edges:
>
> phi.constrain(0., where=mesh.exteriorFaces)
>

> This is actually a Dirichlet (set value) boundary condition. To set
no-gradient (which is no flux in your system), for example, use
phi.faceGrad.constrain(0., where=mesh.exteriorFaces)
The same USAGE page above is helpful here.

>
> I am using the equation:
>
> eq = TransientTerm() == DiffusionTerm(coeff=D) + source - sink
>
>
> I then used the iterative approach as you suggested, and for this setup
> obtain the following after about 60 iterations:
>
> [image: Inline image 1]
>
> Which is a bit surprising. As I only have a source I was expecting the
> whole environment to "flood" with substrate but maybe this is just due to
> me having misunderstood how this PDE is working?
>

First, the flooding (or not) is easier to see if you don't pass
datamax/datamin to the viewer creation. Then, your colobar will change in
time to reflect the _actual_ max and min concentrations you're seeing.
However, in this case, it wouldn't flood indefinitely because the boundary
concentrations are fixed at 0, so you're removing stuff from your system at
the exterior faces. This is easier to see both with a properly scaled
viewer and also if you refine your mesh a bit.
Lx = Ly = 10.
nx = ny = 50
dx = Lx / nx
dy = Ly / ny

[Dario Panada]

Hi Ray,

Thank you, I changed the boundary conditions to

phi.faceGrad.constrain(0., where=mesh.exteriorFaces)

as you suggested, and set my sources and sinks as

source.setValue(50.0, where=(X < 2.0) & (X > 1.0) & (Y > 1.0) & (Y < 2.0))
sink.setValue(20.0, where=(X > 50.0) & (X < 51.0) & (Y > 50.0) & (Y < 51.0))
source.setValue(50.0, where=(X > 95) & (X < 98.0) & (Y > 95) & (Y < 98))

so as to have two sources in bottom left and top right and a sink in
the center. Yes, you are right, what I wait is an infinite source, not
to constraint the value of phi.

With this setup, I'd expect a gradient of substrate to form from the
sources to the center. As the sum of the sources (100.0) is much
greater than the value of the sink, I'd expect substrate to dominate.
But, this is the scenario I get:

[image: Inline image 1]

Where it seems to me like substrate is accumulating on the source and
not diffusing... Any idea why this might be?

Kind Regards,

Dario

[Raymond Smith]

Please continue to include the full script. It appears you've also made
some changes to your original discretization. Nevertheless, I think what's
going on here is a ratio of areas. The source is a "per volume" term. That
is
dc/dt = (diffusion) + (source)
so that the dimensions of each term are something like amount/(volume *
time). So, then, the total amount produced by a given source here (per unit
depth into the screen) is the source value * its area. Note that the source
on the top-right of your system has an area of 16 compared to 1 for the
bottom-left and 1 for the sink in the middle. The result is spreading of
material from the top right corner. When I run a similar simulation, it
does indeed spread from the source. But, with a constant source term, it
will also accumulate in that region simply because diffusive spreading
can't "keep up" with a constant source term.

Ray

[Dario Panada]

Hi,

Apologies, this is the full script http://pastebin.com/eTkFKMX9

I've made a few more changes, and it seems to behave correctly.

I think it might be an issue of visualization because with:

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

[image: Inline image 1]

I can see the substrate spreading and eventually surrounding the sink
(without covering it as the sink consumes more than the sum of the
sources), but with

viewer = Viewer(vars=phi)

I simply see it accumulating?

[image: Inline image 2]

The only thing that's changed is the viewer. I assume what is happening is
that the source in the top right is adding substrate at a speed that
everything else on the surface becomes negligible in comparison and
therefore cannot be seen clearly?

Another thing I noticed - by having sinks, I have some regions which will
have a "negative concentration" until substrate diffuses to them. I was
wondering whether:

a) There was a way to force 0 as a lower limit to the value any grid
position could take;
b) If that made even sense from a "real world" point of view.

In the problem I am trying to model, sinks would be biological cells which
metabolize substrate. If we imagine a metabolic rate per time-step of *m*,
then given a cell with available substrate *a* we have:

*a > m* -> The cell consumes *m* (nutrient requirements fully met)
*0 < a < m* -> The cell does not have enough substrate, so the
concentration at that position becomes 0 and the cell starts starving
*a = 0* -> There is no substrate, the cell starves.

In any case, it should be possible to have a negative concentration. I was
considering sweeping through the array and replacing any negative values
with 0, but I was wondering if there was a more efficient way of going
about this?

Kind Regards,
Dario

[Raymond Smith]

On Sun, Apr 3, 2016 at 3:53 PM, Dario Panada <dario.panada at gmail.com> wrote:

> Hi,
>
> Apologies, this is the full script http://pastebin.com/eTkFKMX9
>
> I've made a few more changes, and it seems to behave correctly.
>
> I think it might be an issue of visualization because with:
>
> viewer = Viewer(vars=phi, datamin=0., datamax=1.)
>
> [image: Inline image 1]
>
> I can see the substrate spreading and eventually surrounding the sink
> (without covering it as the sink consumes more than the sum of the
> sources), but with
>
> viewer = Viewer(vars=phi)
>
> I simply see it accumulating?
>
> [image: Inline image 2]
>
> The only thing that's changed is the viewer. I assume what is happening is
> that the source in the top right is adding substrate at a speed that
> everything else on the surface becomes negligible in comparison and
> therefore cannot be seen clearly?
>

Yes. This.

>
> Another thing I noticed - by having sinks, I have some regions which will
> have a "negative concentration" until substrate diffuses to them. I was
> wondering whether:
>
> a) There was a way to force 0 as a lower limit to the value any grid
> position could take;
> b) If that made even sense from a "real world" point of view.
>

It would probably be better to write your sinks so that they don't consume
if nothing is present. A very simple model which respects this is to say
dc/dt = -k*c
such that the rate of consumption is proportional to the amount already
present at that location. This is referred to as a first-order reaction,
whereas you are modeling the system as a zero-order reaction. The
zero-order approximation _must_ break down as concentration goes to zero,
and there are a number of ways of treating that. You would probably use a
different model than the simple first-order one proposed above, but the
point is that it's probably better to focus on getting a physically
relevant model rather than seeking numerical tricks to bound concentration,
which leads to other problems like not conserving mass within the system.

>
> In the problem I am trying to model, sinks would be biological cells which
> metabolize substrate. If we imagine a metabolic rate per time-step of *m*,
> then given a cell with available substrate *a* we have:
>
> *a > m* -> The cell consumes *m* (nutrient requirements fully met)
> *0 < a < m* -> The cell does not have enough substrate, so the
> concentration at that position becomes 0 and the cell starts starving
> *a = 0* -> There is no substrate, the cell starves.
>
> In any case, it should be possible to have a negative concentration. I was
> considering sweeping through the array and replacing any negative values
> with 0, but I was wondering if there was a more efficient way of going
> about this?
>

You could also potentially have one field variable for nutrients and
another for active cells which can die (but probably not diffuse). Then
they could be coupled reasonably, but this part of the modeling is probably
going to be more up to you and your research.

[Dario Panada]

Hi Ray,

Yes - that makes sense. I think I am now a lot more comfortable with FiPy.

Thank you very much for your help!

Kind Regards,
Dario

[Guyer, Jonathan E. Dr. (Fed)]

Ray -

Thanks for your very complete and very correct answers to Dario.

If you wouldn't mind transcribing your answer to Dario's question on StackOverflow, I'd be happy to up-vote it.

In answer to *your* questions:

> On Apr 3, 2016, at 11:09 AM, Raymond Smith <smithrb at mit.edu> wrote:
>
> Actually, I'm not sure how FiPy treats the steady-state initial guess for Laplace's equation with no flux boundary conditions like yours here. The governing equation + BC's without an initial condition admits any uniform profile as a solution.

> I'm still unsure about the treatment of the initial phi values in this sense as mentioned above.
> 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.

The situation is as you expect. We talk about this toward the end of

http://www.ctcms.nist.gov/fipy/examples/diffusion/generated/examples.diffusion.mesh1D.html

(search for "Fully implicit solutions are not without their pitfalls"). Basically, a steady-state diffusive problem will "lose" its initial condition and should instead be solved by relaxation. The timestep can be made very large; there just needs to be a TransientTerm.

Note, also, that it's not necessary to constrain the gradient to zero to get no-flux. FiPy is a cell-centered finite volume code, so no-flux is the natural boundary condition if nothing else is specified.

- Jon

[Raymond Smith]

Thanks, Jon.

I'll post something to StackOverflow in the next few days.

Also, that makes sense about the steady state ignoring the initial
condition and finding the zero solution. I'd forgotten seeing that in the
example.

Ray

[Raymond Smith]

I posted a reply focusing on the parts of this thread which are in the
original question (and thus on the SO site). Feel free to edit/suggest
edits.

Ray

[Guyer, Jonathan E. Dr. (Fed)]

Looks great, thanks.