Loop macros

[johnomotani]

This PR introduces the loop macros discussed in #35.

This new implementation of parallelization is more flexible - it allows any number of processes to be used, regardless of the number of species or spatial/velocity-space grid dimensions. The dimensions are divided up to give the 'best possible' load balance, with the restriction that each dimension is divided up independently (i.e. the ranges for inner loops are not allowed to depend on the indices of the outer loops).

There is a slight performance drop when evolve_upar == false. The previous implementation had a special case for looping just over ion species, which I have not implemented here to keep things as simple as possible, and assuming that the main target is with evolving upar and ppar anyway. This was used in vpa_advection!(), but the new implementation just tells some processes to continue if they are handling neutrals, when evolve_upar == false.

The implementation in looping.jl is a bit of a nightmare of code-generation, although it's only 375 lines. It was designed to be flexible enough that the only thing we ever have to change is the Tuple of dimensions const all_dimensions = (:s, :z, :vpa) near the top. Sets of loop ranges and macros for loops over any combination of those dimensions are then automatically generated. So hopefully we do not have to fight with it often in future.

Description (copied from README.md)

• There are two types of loop macros:
  • 1. If there is just a single nested loop, then for example
    @s_z_loop is iz begin for izpa in 1:vpa.n f[ivpa,iz,is] = ... end end
    The dimensions in the prefix before _loop give the dimensions that are looped over by the macro, the arguments before begin are the names of the loop variables, in the same order as the dimensions in the prefix; the first dimension/loop-variable corresponds to the outermost nested loop, etc.
  • 2. For more complex loops, a separate macro can be used for each level, for example
    @s_z_loop_s is begin some_setup(is) @s_z_loop_z iz begin @views do_something(f[:,iz,is]) end @s_z_loop_z iz begin @views do_something_else(f[:,iz,is]) end end
    The dimensions in the prefix before _loop_ again give the dimensions that are looped over in the nested loop. The dimension in the suffix after _loop_ indicates which particular dimension the macro loops over. The argument before begin is the name of the loop variables.
• The ranges used are stored in a LoopRanges struct in the Ref variable loop_ranges (which is exported by the looping module). Occasionally it is useful to access the range directly. For example the range looped over by the macro @s_z_loop_s is loop_ranges[].s_z_range_s (same prefix/suffix meanings as the macro).
  • The square brackets [] are needed because loop_ranges is a reference to a LoopRanges object Ref{LoopRanges} (a bit like a pointer) - it allows loop_ranges to be a const variable, so its type is always known at compile time, but the actual LoopRanges can be set/modified at run-time.
• It is also possible to run a block of code in serial (on just the rank-0 member of each block of processes) by wrapping it in a @serial_region macro. This is mostly useful for initialization or file I/O where performance is not critical. For example
  @serial_region begin # Do some initialization f .= 0.0 end
• In any loops with the same prefix (whether type 1 or type 2) the same points belong to each process, so several loops can be executed without synchronizing the different processes. It is (mostly) only when changing the 'type' of loop (i.e. which dimensions it loops over) that synchronization is necessary, or when changing from 'serial region(s)' to parallel loops. To aid clarity and to allow some debugging routines to be added, the synchronization is done with functions labelled with the loop type. For example begin_s_z_region() should be called before @s_z_loop or @s_z_loop_* is called, and after any @serial_region or other type of @*_loop* macro. begin_serial_region() should be called before @serial_region.
• Internally, the begin_*_region() functions call _block_synchronize(), which calls MPI.Barrier(). When all debugging is disabled they are equivalent to MPI.Barrier(). Having different functions allow extra consistency checks to be done when debugging is enabled, see debug_test/README.md.

Closes #35, closes #38.