Lagrangian particle tracking using eval

[sijiehuang23]

Hi,

I see in the source code utilities/lagrangian_particle.py that the shenfun code can interpolate fields using the eval method from shenfun.Function.

However, the example is rather small (a $40^2$ flow field). So I'm wondering, does this work in an MPI context (when the simulation is running in parallel)? Is it efficient? Can I use this functionality for large-scale particle tracking (for example, in a $1024^3$ simulation)?

[mikaem]

Yes, the particle approach works in parallel. At least it should, haven't tried it in a while.

[sijiehuang23]

I wrote a code to test it with 4 MPI processes. And yes, I think you're right—it does support parallelization with MPI (at least on the face of it).

I'm showing the result of a 2D random field (here I only consider interpolation points at constant $y$). The result from shenfun is quite good, compared with a cubic spline from scipy.

However, eval requires the full data array, instead of distributed ones. In other words, each process has a copy of the full input data array. It does not work with input position array x distributed among processes.

[sijiehuang23]

On top of this, I found thie evaluation becomes incredibly slow for large number of particles (say $\gtrsim10^5$). And increasing number of processes won't help.

[mikaem]

You are right @sijiehuang23, the particle code is not very efficient. One of the reasons for this is that you cannot use the FFT for the particles, because the FFT only works for uniform meshes. As of now the particles velocity is computed with full spectral accuracy using eval and for this reason the particles position must be available on all processes. This is ok for a relatively small number of particles, but it does not scale well, as pointed out. However, note that since eval is computed in spectral space it does not really make sense to let the particles live only on the process where their position is at. This makes a lot of sense for a low order finite difference/volume/element code, but not for fully spectral. I could envision, though, a code where a particle's velocity is computed locally to lower order using only neighbouring points in physical space and interpolation.

[sijiehuang23]

Yes, I think you're right. As far as I know, most particle codes use polynomial interpolation in physical space. Though some uses non-uniform FFT.