perrette
commented
9 years ago After a bit of experimenting:
- fortran + f2py seem by far the easiest option to generate fast code !!
- https://github.com/pism/regional-tools is a really handy C++ library to label basin upslope a few pixels. The underlying code is pretty simple, but one key step of building a streamline based on a vector field (in that case, downward elevation gradient) is done with the gsl_odeiv C++ library using Runge-Kutta-Fehlberg algorithm, instead of simply adding a fixed-step vector to current position. This is supposed to minimize numerical errors when integrating the vector field, and indeed seems to work much better. An option could be to re-use bits of code from pism / regional-tools to derive the streamline in a similar way.
Cost associated with adopting C++ code from pism-regional tool:
- more complex C++ (2-D array is not a native type, use of pointers)
- cython wrapper quite verbose, instead of just running f2py...
- the above two points may not be critical because the basic pism-regional tool already exists, so it is relatively easy to marginally modify the code. It may become problematic at some point if further extension of the code are required, and it makes it more difficult (up to impossible) to reuse code into existing fortran code.
Issues / advantages with using a fortran implementation
- possible need to reprogram the Runge-Kutta.Fehlberg algorithm, or to understand existing libraries
- reuse in existing fortran code
Speed?
- I have not performed speedtest between a fortran and a C++ implementation yet
- within the C++ code, the Runge-Kutta-Fehlberg method yields 50% increase in execution time compared to a more naive approach.
Flowline calculation is currently relatively slow. One part is related to necessary loading of relatively large netCDF files, the other relates to the calculation itself. The first - because easier - task would be to speed-up the calculation part, for instance by writing a very short fortran code with iso_c_bindings (or using f2py), or with cython, or whatever. It is unclear to me whether the loading of netCDF files could be made quicker, since netCDF4 is already a wrapper around netCDF4/HDF5 routines.
An additional option could be to propose the drawing of a series of flowlines which could be seeded from an existing polyline, so that the netCDF file is loaded one, and all flowlines calculated at once.