seismic-anisotropy / PyDRex

Simulate crystallographic preferred orientation evolution in polycrystals
https://seismic-anisotropy.github.io/PyDRex/
GNU General Public License v3.0
3 stars 2 forks source link

Add Fortran routines for performance comparisons #208

Open adigitoleo opened 1 month ago

adigitoleo commented 1 month ago

Maybe the best way to do it is using numpy.f2py to generate some Fortran extensions. This way, we can pass data from Python (much nicer for I/O and setup) and even run the Fortran stuff in CI as part of the test suite if desired. Manuele Faccenda has a refactored version of DRex with a more generic interface compared to the original. It seems like the newest Fortran DRex variant and is actively used, which makes it the best candidate for comparisons. I think I'll use that code as the basis for our Fortran extensions.

It seems like the preferred way of writing Fortran extensions that gives the most control over the data transfer involves adding f2py directives in the Fortran source (just as writing C extensions actually involves writing "Cython").

We want drop-in replacements for the following PyDRex methods:

We also might want the statistical methods:

And it could be nice to have:

adigitoleo commented 1 month ago

Relevant:

Maybe better than f2py:

Lfortran is working on their own thing as well

Distributing any of this with the PyPI wheel is a non-goal, too many hoops to jump through. Only needs to work on Linux (gfortran, probably).

adigitoleo commented 1 month ago

Starting work on the Fortran routines in add-fortran-routines (branch). Aiming for Fortran 2003 compatibility which provides proper inf values and real64, among other things. gfortran -fsyntax-only is a good way to lint, as well as fortls which I added to the dev dependencies.

To call the functions/subroutines from Python, fmodpy looks the most appealing of the above options, but I have yet to test it properly (passing in numpy arrays, GC, parallel? or we only test serial performance?). f2py probably only supports up to Fortran 95, and going through a C layer is the more modern approach. fmodpy also doesn't require any special annotations in the Fortran source.