speeding up near2far

[stevengj]

Two options:

• use #pragma omp parallel for in farfield_lowlevel, maybe by moving for (int i = 0; i < Nfreq; ++i) to be the outermost loop. This only works if you have cores to spare and we compile with OpenMP. openmp for near2far calculation #868

• in 2d, we can switch to a faster Bessel-function calculation. In a quick test, the Fortran routines by Amos libc jn and yn routines are about 4-5x faster than the ones in GSL. Precompiled binaries of these are available in https://github.com/JuliaMath/openspecfun at https://github.com/JuliaMath/OpenspecfunBuilder/releases, for example, if we don't want to require a Fortran compiler. move away from GSL functions. #869

For the first option, would be nice to try it out. Note that you will need to compile with -fopenmp

[stevengj]

I added a branch omp_near2far with the changed loop (871b4ca). You will need to manually add the -fopenmp flag by configuring with

./configure CXX="g++ -fopenmp"

or similar.

Then you set the number of threads with the environment variable OMP_NUM_THREADS. For example:

OMP_NUM_THREADS=4 meep somefile.ctl

[oskooi]

This seems to be working and giving rise to linear speedup as demonstrated in the figure below. The test involves timing the get_farfields calculations in python/examples/binary_grating_n2f.py involving 21 frequencies and 500 points for 5 values of OMP_NUM_THREADS (1-5) for two cases: (1) a unit cell with add_near2far involving nperiods=20 and (2) a super cell with 41 unit cells and nperiods=1 (i.e., default: no tiling). Note: since the total time for get_farfields is not yet directly available, pending #856, the values were computed by subtracting from the "Elapsed run time" all the components in "Field time usage" (e.g., "time stepping, "communicating", etc).

The only unusual thing was that it worked for the serial version (i.e., CXX="g++ -fopenmp") but not with the OpenMPI compiler (i.e., CXX="mpic++ -fopenmp").

[stevengj]

If you use MPI, you need to set MPICXX="mpic++ -fopenmp" because of this line

[stevengj]

Actually, the Fortran Amos routines are just for complex arguments, which we don't need. For real arguments and integer order, Julia is using the POSIX bessel functions, which we should use too.