exafmm Notes

exafmm-notes.txt

Notes on using exafmm for multiple evaluation on moving particles and common pitfalls.

1) NCRIT
You need to choose NCRIT >= 2 even for debugging, because else too little memory is allocated.

2) Store Bodies
The Cell structs use pointers into the Bodies vector. This means the Bodies vector needs to be destructed after Cells,
i.e. as global variable or class member.

3) Ordering
exafmm reorders Bodies during calculation, this means for tracing particles, we need to add an index. This can be implemented like this:

struct Body {
  exafmm::vec3 X;                                     //!< Position
  exafmm::real_t q;                                   //!< Charge
  exafmm::real_t p;                                   //!< Potential
  exafmm::vec3 F;                                     //!< Force
  size_t index;
  bool operator<(const MyBody &other) { return index < other.index; };
};

Then we can use std::sort(bodies.begin(), bodies.end()) to sort them afterwards.

4) Multiple evaluation
exafmm supports multiple evaluation of the field of the initial Bodies on a set of test Bodies.
Given: charged_bodies, test_bodies

// global variables
NCRIT = 64;
THETA = 0.4;
P = 10;

// initialize charged Bodies
initKernel();
Cells charged_cells = buildTree(charged_bodies)
upwardPass(charged_cells);

for(uint step=0; step < num_iterations; step++) {
    Cells test_cells = buildTree(test_bodies);
    for(size_t i=0; i < test_cells.size(); i++) {
        test_cells[i].L.resize(NTERM);
        test_cells[i].M.resize(NTERM);
    }
    horizontalPass(test_cells, charged_cells);
    downwardPass(test_cells);
    std::sort(test_cells.begin(), test_cells.end());

    // Do something, i.e. moving test_bodies[i] into test_bodies[i].F direction
}

5) Array definition with variable size
Some compilers do not support C array allocation with non-constant size in kernel.h:
complex_t Ynm[P*P], YnmTheta[P*P];

Create global static variables in kernel.h:
static std::vector<complex_t> Ynm, Ynm2, YnmTheta;

add to initKernel():
Ynm.resize(P*P);
YnmTheta.resize(P*P);
Ynm2.resize(4*P*P);

6) OpenMP 3.0
Older compilers and MSVC do not support OpenMP 3.0. Comment out all
#pragma omp task
#pragma omp taskwait
directives in traverse_eager/lazy.h. The
#pragma omp parallel
#pragma omp single nowait
directives are okay.

7. buildTree will crash, if bodies.size() is 0.