Julia调用PyFR

PyFR由Python3写成，并且不提供接口，所以我们需要自己动手，改造一番，才能在Julia里调用它的内部函数。上回我们介绍了PyFR的运行示例：二维欧拉涡。在该示例中，有一个转换.msh文件的步骤。如何在Julia调用PyFR实现这个步骤呢？我们将从这里开始，最终实现Julia调用PyFR运行完整仿真。一切开始前，假定你已经安装了Julia的PyCall包，并且读过PyCall的Github项目的readme文档。

入门：如何调用PyFR函数转换网格文件

首先，我们要安装PyFR。这次不是在Python虚拟环境里，而是正常环境里安装，直接在终端执行：sudo apt install pyfr。

安装后，仍然需要上回提到的源代码中的examples，切换到其中的euler_vortex_2d目录下。尝试运行：pyfr import euler_vortex_2d.msh euler_vortex_2d.pyfrm，看看是否正常运行。

被安装的PyFR的源代码在/usr/local/lib/python3.7/dist-packages/pyfr/里，打开__main__.py文件，里面有几个至关重要的函数。_process_common是最重要的求解函数，但在了解它之前，需要先学习怎么调用一些简单的函数，比如我们现在要讲的process_import，它接收一个参数args，并使用args中的变量执行若干操作，将.msh文件转换为.pyfrm。

args实际上是一个Namespace，是python中的一种特殊的字典。我们注意到，__main__.py里的函数多以args为输入参数，所以很有必要搞清楚args到底包含了什么内容。文件开头有个main函数，利用python的argparse包来自定义终端命令。仔细看看，这段代码也是学习argparse的优秀样例。其中有一行args = ap.parse_args()，正是创建args的地方。我们在下方添加一行：

    args = ap.parse_args()
    print(args)

并保存文件（需要权限）。

然后在终端里重新运行pyfr import euler_vortex_2d.msh euler_vortex_2d.pyfrm，会看到如下输出：

Namespace(cmd='import', inmesh=<_io.TextIOWrapper name='euler_vortex_2d.msh' mode='r' encoding='UTF-8'>, outmesh='euler_vortex_2d.pyfrm', process=<function process_import at 0x7f4ab21be7a0>, type=None, verbose=None)

这就是args的真面目。当我们运行PyFR的其他命令如pyfr run时，会看到输出的这个Namespace里包含的内容是不同的。

Julia里貌似没有对应的Namespace。用PyCall虽然可以调用argparse包，但无法使用自定义命令。所以我们应当舍弃Namespace，把args分解为多个变量，传给函数以达成目的。

具体来说，我们观察process_import里的语句，发现此处的args有三个变量：type, inmesh, outmesh。其中type可以为[]，inmesh是python读取.msh文件创建的io流，outmesh是转换文件的文件名。于是我们把函数重新写成：

def process_import(type,inmesh,outmesh):
    # Get a suitable mesh reader instance
    if type:
        reader = get_reader_by_name(type, inmesh)
    else:
        extn = os.path.splitext(inmesh.name)[1]
        reader = get_reader_by_extn(extn, inmesh)

    # Get the mesh in the PyFR format
    mesh = reader.to_pyfrm()

    # Save to disk
    with h5py.File(outmesh, 'w') as f:
        for k, v in mesh.items():
            f[k] = v

保存。打开另一个文件__init__.py。添加两行：

from pyfr.__main__ import process_import as pimp
__all__ = ['pimp']

这个__init__.py文件就是python package的导出接口，凡是出现这个作用域里的东西都可以被外部引用。我们添加的第一行把__main__文件里的process_import函数导入该作用域并命名为pimp，这样便可以在外部引用它了。__all__变量只是一个提示用的名称列表，可以不写。保存。

现在我们在Julia REPL里操作。首先在euler_vortex_2d目录里打开新的终端并进入Julia REPL。

第一步是调用必要的包。

using PyCall
pyfr = pyimport("pyfr")

我们用python的io包来创建inmesh变量：

io = pyimport("io") 
inmesh = io.open("euler_vortex_2d.msh","r",encoding="utf-8")

声明一个type和outmesh变量：

type = []
outmesh="euler_vortex_2d.pyfrm"

调用pimp函数（注意文件路径）：

pyfr.pimp(type,inmesh,outmesh)

可以看到文件夹里新生成的.pyfrm文件，表示已完成。

最后我们可以查看一下刚才写的那个列表，看看有哪些东西可调用：

julia> pyfr.__all__
1-element Array{String,1}:
 "pimp"

小贴士一：注意__main__.py 文件前后都是两条下划线。有时下划线是一条，有时是两条，容易看错。
小贴士二：每次修改pyfr源文件后，不会自动在python环境或Julia REPL环境里更新，需要重新import，实在不行就重启环境。

进阶：如何调用PyFR函数完成仿真

有了以上基础后，剩下的皆水到渠成。我们在euler_vortex_2d目录里打开终端，运行：

pyfr run -b openmp -p euler_vortex_2d.pyfrm euler_vortex_2d.ini

之前我们在__main__.py里写的那句print(args)此刻再立新功，为我们揭示了pyfr run命令下的args内容：

Namespace(backend='openmp', cfg=<_io.TextIOWrapper name='euler_vortex_2d.ini' mode='r' encoding='UTF-8'>, cmd='run', mesh='euler_vortex_2d.pyfrm', process=<function process_run at 0x7fafd342f680>, progress=True, verbose=None)

其中对我们有用的是backend、cfg、mesh三个参数。backend和mesh都是字符串，cfg是io流。延续之前的思路，我们把args分解，将实现pyfr run命令的两个函数_process_common和process_run改写为：

def _process_common(backend, progress, mesh, soln, cfg):
    # Prefork to allow us to exec processes after MPI is initialised
    if hasattr(os, 'fork'):
        from pytools.prefork import enable_prefork

        enable_prefork()

    # Import but do not initialise MPI
    from mpi4py import MPI

    # Manually initialise MPI
    MPI.Init()

    # Ensure MPI is suitably cleaned up
    register_finalize_handler()

    # Create a backend
    backend = get_backend(backend, cfg)

    # Get the mapping from physical ranks to MPI ranks
    rallocs = get_rank_allocation(mesh, cfg)

    # Construct the solver
    solver = get_solver(backend, rallocs, mesh, soln, cfg)

    # If we are running interactively then create a progress bar
    if progress and MPI.COMM_WORLD.rank == 0:
        pb = ProgressBar(solver.tstart, solver.tcurr, solver.tend)

        # Register a callback to update the bar after each step
        callb = lambda intg: pb.advance_to(intg.tcurr)
        solver.completed_step_handlers.append(callb)

    # Execute!
    solver.run()

    # Finalise MPI
    MPI.Finalize()


def process_run(backend, progress, mesh, cfg):
    _process_common(
        backend, progress, NativeReader(mesh), None, Inifile.load(cfg)
    )

简单地说，就是把参数改一改。注意到这里有个progress参数，一般设为true即可，用于显示一个仿真的进度条。

我们在外部调用的是process_run，所以在__init__.py里写成：

# -*- coding: utf-8 -*-

from pyfr._version import __version__
from pyfr.__main__ import process_import as pimp
from pyfr.__main__ import process_run as prun

__all__ = ['pimp','prun']

最后，我们在Julia里调用prun来执行，不妨写个脚本：

using PyCall
io = pyimport("io")
pyfr=pyimport("pyfr")
mesh = "euler_vortex_2d.pyfrm"
cfg = io.open("euler_vortex_2d.ini","r",encoding="utf-8")
# 运行仿真
pyfr.prun("openmp",true,mesh,cfg)

应该能看到进度条。恭喜你完成了。