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?
  • $\begingroup$ What you describe reminds me (somewhat) of Particle in Cell (PIC) Methods: en.wikipedia.org/wiki/Particle-in-cell $\endgroup$
    – Paul
    Sep 13, 2014 at 14:38
  • $\begingroup$ Fenics is a finite element package. $\endgroup$
    – Jan
    Sep 13, 2014 at 16:56
  • $\begingroup$ 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$? $\endgroup$ Sep 14, 2014 at 6:09
  • $\begingroup$ Look up "particle-mesh methods". $\endgroup$ Sep 14, 2014 at 6:12
  • $\begingroup$ @Paul yes it is essentially PIC, just thought that I'd minimise jargon. $\endgroup$ Sep 14, 2014 at 14:50

1 Answer 1


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

In general, it is not entirely trivial to find out which cell a particular point lies in if you have an unstructured mesh. There are, however, a number of simpler cases:

  • If you have have a mesh composed of triangles (in 2d) or tetrahedra (3d), then each cell is bounded by linear constraints and finding whether a point is inside a cell only requires to check 3 (in 2d) or 4 (in 3d) linear inequalities.

  • If you have quadrilateral meshes in 2d, each cell is bounded by 4 linear inequalities as long as every cell is convex (which you need for other reasons as well).

  • In 3d with hexahedra, cells are bounded by nonlinear constraints and things become a lot more complicated.

But even in the case of linear constraints, one surprisingly frequently runs into the case where points are within round-off distance from a bounding line/plane and the function that determines whether the point lies in a cell may say "yes" for two or no cell. You will have to figure out how to deal with this.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.