Calling Fortran from Python using CFFI: using WRF code from Python

Why do you want to use Fortran code?

Fortran is already a forgotten language by most software developers, but it is still widely used in some fields of science; atmospheric science being a good example. Therefore, you have a large amount of software already written and it’s still growing…
If you work on really big data and if the performance is (almost) the only priority you have. Fortran is still the fastest language for some applications, although the advantage over other languages is often overestimated, but that’s another story.

Why do you want to call from Python?

It’s easier for you to write additional code using Python than Fortran.
You would like to use Python environment, e.g., for testing purposes.

CFFI – C Foreign Function Interface for Python

Enables you to call compiled C and Fortran code from Python.
Enables you to call Python code from C and Fortran using callback mechanism.
Attempts to support both PyPy and CPython.

How to do it – step by step

Creating a shared library

You have to first create a shared library from your Fortran code.

You can do it by compiling your code with -shared option, e.g.:

gfortran -shared  fortran_code.f90 -o libfortran.so  -fPIC

I understand that different fortran compilers can create different function’s signatures in the library and that they are not easy to guess, so you should check the names using the nm command. In addition, OSX can add one additional underscore you should skip later.

You can also use iso_c_binding module to write a wrapper to your Fortran function or subroutine. This approach requires additional coding, but:
- forces Fortran compiler to produce the correct C signature
- allows passing arguments as values (not by references)
- allows to insert code that will handle conversion to Fortran-specific argument types (e.g. deferred-length strings)
- helps to assure proper argument types

Creating a wrapper

In the example presented in the post, writing a wrapper was probably not absolutely needed, but I decided to write it anyway, so will explain briefly how to do it.
Your fortran_wrapper.f90 should import the fortran module you want to wrap and all c-types from iso_c_bindings you will use

use your_fortran_module
use iso_c_binding, only: c_int, c_double

You have to define a C fuction/subroutine – binding with the Fortran function/subroutine you want to use:

 subroutine c_wrapper_your_subr(arg1, arg2, arg3) bind(c)

You should declare all variables that will be passed as values – arg1; and as pointers – arg2, arg3:

integer(c_int), intent(in), value :: arg1  
real(c_double), intent(inout) :: arg2, arg3(1:10)

Now, you can call the Fortran function/subroutine from your_fortran_module:

call your_subr(arg1, arg2, arg3)

Having a wrapper to your Fortran module, you can create a shared library:

gfortran -c -fPIC fortran_code.f90 -o fortran_code.o
gfortran -shared -fPIC  fortran_code.o fortran_wrapper.f90 -o libfortran.so

Calling Fortran function/subroutine from Python

Once you have a shared library, you can call your subroutine from Python. First, you need to load the library using a foreign function interpreter FFI:

from cffi import FFI
ffi = FFI()
lib = ffi.dlopen("libfortran.so")

For each function/subroutine you want to call from Python, you should provide a C-signature to CFFI:

ffi.cdef("void c_wrapper_your_subr(int arg1, double arg2, double arg3[]);")

Before calling your subroutine you should define arguments. If you’re passing just values, you can simply write:

arg1 = 1
arg2 = 3.14

Since arg3 is an array, you should create a CFFI object storing pointer to the NumPy array. Note that NumPy array should have a Fortan order (however for 1d it doesn’t change anything) :

numpy_array = numpy.arange(10, order='F')
arg3  = ffi.cast("double*", numpy_array.__array_interface__['data'][0])

Now, you can call your Fortran subroutine:

lib.c_wrapper_your_subr(arg1, arg2, arg3)

Example of calling a WRF microphysical scheme

In my work I wanted to test a new C++ library (libcloudphxx.igf.fuw.edu.pl ) of microphysical schemes. As a part of testing, I decided to compare results for simple cases with results obtained with a widely used microphysical schemes from the Weather Research and Forecasting (WRF; http://www.wrf-model.org/index.php) model that are written in Fortran. I decided to call various microphysical schemes from WRF and the new C++ library using Python and compare them for exactly the same setup.

Below I’m presenting an example of calling the simplest WRF microphysical scheme based on the Kessler approach. All presented codes are on github repository - https://github.com/djarecka/scientific-software-diary. The WRF microphysical scheme, module_mp_kessler.f90, is taken from http://www2.mmm.ucar.edu/wrf/users/download/get_sources.html

The wrapper for this case, is quite long, because the Fortran subroutine has a long list of arguments. Some of arguments, I’m passing as values. All arrays will be passed as pointers, most of them (temperature, water vapor and cloud/rain mixing ratios) will be changed by the Fortran subroutine. The full content of kessler_wrap.f90:

module kessler_wrap
  use iso_c_binding, only: c_int, c_double
  use module_mp_kessler
  implicit none
contains
 
  subroutine c_kessler(t, qv, qc, qr, rho, pii ,dt_in, z, xlv, cp, &
                       EP2, SVP1, SVP2, SVP3, SVPT0, rhowater,     &
                       dz8w, RAINNC, RAINNCV,                      &
                       ids,ide, jds,jde, kds,kde,                  &
                       ims,ime, jms,jme, kms,kme,                  &
                       its,ite, jts,jte, kts,kte                   &
                      ) bind(c)
 
    real(c_double), intent(in), value :: dt_in,  xlv, cp, rhowater,    &
                                         EP2,SVP1,SVP2,SVP3,SVPT0
    integer(c_int), intent(in), value :: ids,ide, jds,jde, kds,kde,    &
                                         ims,ime, jms,jme, kms,kme,    &
                                         its,ite, jts,jte, kts,kte
    real(c_double),  intent(in)       :: rho(ims:ime,kms:kme,jms:jme), &
                                         pii(ims:ime,kms:kme,jms:jme), &
                                         dz8w(ims:ime,kms:kme,jms:jme),&
                                         z(ims:ime,kms:kme,jms:jme)
    real(c_double),  intent(inout)    :: t(ims:ime,kms:kme,jms:jme),   &
                                         qv(ims:ime,kms:kme,jms:jme),  &
                                         qc(ims:ime,kms:kme,jms:jme),  &
                                         qr(ims:ime,kms:kme,jms:jme),  &
                                         RAINNC(ims:ime,jms:jme),      &
                                         RAINNCV(ims:ime,jms:jme)
 
    call kessler(t, qv, qc, qr, rho, pii, dt_in, z, xlv, cp,   &
                 EP2, SVP1, SVP2, SVP3, SVPT0, rhowater,       &
                 dz8w, RAINNC, RAINNCV,                        &
                 ids,ide, jds,jde, kds,kde,                    &
                 ims,ime, jms,jme, kms,kme,                    &
                 its,ite, jts,jte, kts,kte)
 
    end subroutine c_kessler
end module kessler_wrap

Creating a shared library and forcing compiler to use double precision:

gfortran -fdefault-real-8 -c -fPIC module_mp_kessler.f90 -o module_mp_kessler.o
gfortran -fdefault-real-8 -shared -fPIC module_mp_kessler.o kessler_wrap.f90 -o libkessler.so

In the cffi_kessler.py I defined a python function that takes dictionary of all needed arrays with atmospheric variables (variable_nparr), and arrays dimensions as arguments. Since, I have to pass many arrays to Fortran, I have to create dictionary of CFFI objects storing appropriate pointers (variable_CFFI). Variables ims, ime… that are used to allocate enough memory for arrays in the Fortran subroutine I’m passing just as values, so no need to create any additional CFFI objects.

import numpy as np
from constants_kessler import xlv, cp, EP2, SVP1, SVP2, SVP3, SVPT0, rhowater
from cffi import FFI
ffi = FFI()
 
# function creates cdata variables of a type "double *" from a numpy array             
# additionally checks if the array is contiguous                                       
def as_pointer(numpy_array):
    assert numpy_array.flags['F_CONTIGUOUS'], \
        "array is not contiguous in memory (Fortran order)"
    return ffi.cast("double*", numpy_array.__array_interface__['data'][0])
 
def kessler(nx, ny, nz, dt_in, variable_nparr):
    ffi.cdef("void c_kessler(double t[], double qv[], double qc[], double qr[], double rh\
o[], double pii[], double dt_in, double z[], double xlv, double cp, double EP2, double SV\
P1, double SVP2, double SVP3, double SVPT0, double rhowater, double dz8w[], double RAINNC\
[], double RAINNCV[], int ids, int ide, int jds, int jde, int kds, int kde, int ims, int \
ime, int jms, int jme, int kms, int kme, int its, int ite, int jts, int jte, int kts, int\
 kt);", override=True)
 
    # load a library with the C function                                          
    lib = ffi.dlopen('libkessler.so')
 
    # create cdata variables for each numpy array            
    variable_CFFI = {}
    for item in ["t", "qv", "qc", "qr", "rho", "pii", "z", "dz8w", "RAINNC", "RAI\
NNCV"]:
        variable_CFFI[item] = as_pointer(variable_nparr[item])
 
    [ims, ime, ids, ide, its, ite] = [1, nx] * 3
    [jms, jme, jds, jde, jts, jte] = [1, ny] * 3
    [kms, kme, kds, kde, kts, kte] = [1, nz] * 3
 
    # call the C function                                                                        
    lib.c_kessler(variable_CFFI["t"], variable_CFFI["qv"], variable_CFFI["qc"],
                  variable_CFFI["qr"], variable_CFFI["rho"], variable_CFFI["pii"],
                  dt_in, variable_CFFI["z"], xlv, cp, EP2, SVP1, SVP2, SVP3, SVPT0,
                  rhowater, variable_CFFI["dz8w"], variable_CFFI["RAINNC"],
                  variable_CFFI["RAINNCV"], ids, ide, jds, jde, kds, kde,
                  ims, ime, jms, jme, kms, kme, its, ite, jts, jte, kts, kte)

The simplest example of use for the python function is in kessler_call.py that you can find in the repository https://github.com/djarecka/scientific-software-diary with all other codes.

Thanks for help to:

Sylwester Arabas
Davide del Vento
Maciej Fijalkowski