# Numerically solve a PDE in Python with a term calculated by coarse-graining

I'm trying to solve a PDE in Python of the form,

$\dfrac{\partial c(\mathbf{x}, t)}{\partial t} = \mathrm{D} \nabla^2 c(\mathbf{x}, t) -\gamma \rho(\mathbf{x}, t) c(\mathbf{x}, t)$

where $c$ represents the concentration of a chemical, $\mathrm{D}$ and $\gamma$ are constants, and $\rho$ represents the density field of some large point-particles (which in fact model bacteria). This is all done in 2D.

The problem is that I'm modelling the particles/bacteria represented by $\rho$ as individual agents moving in continuous space, not as a density field obeying a PDE.

How I get round this is to take the positions of all the particles and turn this into a coarse density field $\rho$, then use this to solve the PDE.

Up until now I've been doing this by hand, binning particle positions onto a cartesian grid, then iterating the PDE using finite differencing. So far this has worked fine.

But now I want to solve the PDE in a more complicated 2D geometry with curved boundaries, so I think I need to use the finite volume method rather than finite differencing. This is beyond what I can do by hand, so I've been looking into packages like FiPy and FEniCS.

For me, there seem to be two issues I don't know if they can deal with:

1. Is there an easy way to determine in which volume element a particular point lies? If so I can calculate $\rho$ by hand at each time step so that's good.
2. Assuming I can calculate $\rho$, can I supply a variable field like this to a PDE solver, which it will then use in iterating the equation?
• What you describe reminds me (somewhat) of Particle in Cell (PIC) Methods: en.wikipedia.org/wiki/Particle-in-cell – Paul Sep 13 '14 at 14:38
• Fenics is a finite element package. – Jan Sep 13 '14 at 16:56
• So you can compute $\rho$ for all time before starting to solve the evolution equation for $c$? That is, there is no feedback from $c$ into the evolution of $\rho$? – David Ketcheson Sep 14 '14 at 6:09
• Look up "particle-mesh methods". – David Ketcheson Sep 14 '14 at 6:12
• @Paul yes it is essentially PIC, just thought that I'd minimise jargon. – EddieJessup Sep 14 '14 at 14:50

The way it sounds to me is that the key step is really to determine the field $\rho$ by evaluating how many particles you have in each cell. The rest is then just a standard heat equation with a particular source term, which is not overly complicated. For example, you could start at step-26 of the deal.II tutorial: http://www.dealii.org/developer/doxygen/deal.II/step_26.html (disclaimer: I wrote that program and a significant part of the library behind it).