Add Fortran routines for performance comparisons

[adigitoleo]

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:

• minerals.update_all, which requires
• minerals.update_orientations, which requires
• core.derivatives

We also might want the statistical methods:

• minerals.voigt_averages
• stats.resample_orientations (was this ever included in the Fortran?)

And it could be nice to have:

• diagnostics.elasticity_components (but check the new routine from Tape & Tape 2024 first)
• diagnostics.finite_strain

[adigitoleo]

Relevant:

• https://fortran-lang.discourse.group/t/interfacing-fortran-code-from-python/1280/16
• https://www.matecdev.com/posts/fortran-in-python.html

Maybe better than f2py:

• https://github.com/jameskermode/f90wrap (produces wrappers for f2py)
• https://github.com/tchlux/fmodpy (generates a Fortran<->C<->Python pipeline via C wrappers on either side)

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]

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.

[adigitoleo]

Although I was able to get drop-in replacements for the methods in core to run, they took a long time and produced some strange outputs which I haven't had time to interrogate. Note that data transfer with/from the fortran wrappers requires numpy.asfortranarray on the sending side and np.ascontiguousarray on the receiving side (I think). I ran into some bugs in fmodpy itself where it tries to do np.zeros to initialise array sizes but uses a c_int without appending .value to get the normal python integer. Looking at the fmodpy source did not make me excited to contribute fixes there either. It may be worth rolling the C interface by hand if there is interest in making these drop-in routines work.

Due to the difficulties, I am dropping this from the beta release milestone.

[Patol75]

That is a sensible choice from my perspective. In terms of performance comparison, I would be curious to see the difference between PyDRex and the initial implementation. Perhaps we could just use the test case we had years ago; I still have the files on Gadi. What do you think?

[adigitoleo]

The one I used to test for the paper revisions was tools/drex_forward_shear.f90 which skips the pathline computation and just uses deriv and decsym (I compared it to an equivalent 2D simple shear case in Python, which included a warm-up compilation for Numba that was not timed -- I should have a file for that somewhere).