使用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数组指针

Note

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从Python调用Fortran程序:Python调用WRF代码

Contents

为什么要使用Fortran代码

为什么从Python调用Fortran

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

操作步骤

创建共享库

函数封装器

从Python调用Fortran函数/子程序

调用WRF微物理方案

延伸阅读

相关文章