使用CFFI从Python调用Fortran程序:Python调用WRF代码

为什么要使用Fortran代码

• Fortran语言是已经几乎被遗忘的语言，但其在科学领域却仍为活跃，尤其是大气科学领域。
• 如果效率优先是需要考虑的重要因素，那么可以尝试使用Fortran，虽然Fortran在这方面的作用相较于其他语言经常被高估(=_=)。

为什么从Python调用Fortran

• 使用Python比使用Fortran写模块更简单
• 测试方面，使用Python环境更友好

CFFI：面向Python的C语言外部函数接口

• 从Python调用已编译的C和Fortran代码
• 使用回调机制从C和Fortran中调用Python程序
• 试图支持PyPy和CPython

操作步骤

创建共享库

首先创建Fortran程序的共享库

• 可以使用-shared选项编译你的代码

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

  不同的fortran编译器在库中可能会创建不同的函数签名，因此在使用之前应该使用nm命令检查名称。

• 也可以使用iso_c_binding模块对Fortran函数和子程序进行封装，此方法需要额外的代码，但这会：

  • 强制Fortran编译器生成正确的C签名
  • 允许传递参数作为值
  • 允许嵌入代码处理Fortran特定参数类型的转换
  • 有助于确保正确的参数类型

函数封装器

在这个示例中，我们会创建一个封装器，虽然可能并不需要这样。

fortran_rapper.f90应该导入想要封装的Fortran模块，并且从iso_c_binding中调用需要使用的c_types：

1
2

use your_fortran_module
use iso_c_binding, only: c_int, c_double

你需要定义一个C函数/子程序—打包需要使用的Fortran函数/子程序：

1

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

然后需要声明所有需要使用的变量，arg1是值，arg2和arg3是指针：

1
2

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

现在就可以从Fortran模块your_subr中调用上述函数/子程序了：

1

call your_subr(arg1, arg2, arg3)

最后一步，创建共享库封装Fortran模块，

1
2

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

从Python调用Fortran函数/子程序

有了共享库之后，就可以从Python中调用子程序了。

首先，使用外部函数解释器FFI加载库

1
2
3

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

对于那些需要从Python调用的函数/子程序，需要提供C签名给CFFI：

1

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

调用子程序之前需要定义参数，如果只是传递值，可以简单的赋值：

1
2

arg1 = 1
arg2 = 3.14

因为arg3是一个数组，你应该创建一个CFFI对象存储Numpy数组指针

NumPy数组的存储方式应该和Fortran数组相同。

1
2

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

现在就可以调用Fortran子程序了

1

lib.c_wrapper_your_subr(arg1, arg2, arg3)

调用WRF微物理方案

为了测试C++版的微物理方案，需要使用WRF中广泛使用的微物理方案来进行对比。

下面以WRF中Kessler方法为例，所有代码在这里，WRF的微物理方案module_mp_kessler.f90可以从这里下载。

• 因为WRF的微物理方案需要传递大量的参数，封装起来会比较长。所有的数组将直接传递指针，以下是kessler_wrap.f90的内容：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

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                                     

• 强制编译器使用双精度创建共享库

  1 2	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

• 在cffi_kessler.py中，定义一个python函数处理所有需要的大气变量数组，以及数组维度。因为需要传递大量的数组给Fortran，因此必须要创建一个CFFI字典来存储合适的指针。变量ims，ime…用来在Fortran子程序中分配足够的内存给数组：

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40

  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)

• 这是使用Python函数的一个最简单示例，在这里你可以发现更多的代码。

延伸阅读

• CFFI documentation
• CFFI example I used
• WRF model
• Libcloudph++ library
• Fortran signatures
• All examples are on github

原文链接