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

Question

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?

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.