SPHYSICS流体力学仿真模拟程序的动态链接库编译及C#混合编程方法

想写的原因很简单，怕时间长了，自己也不记得过程和方法了。正好可以写下来，以后和有这方面需要的朋友们探讨。

首先说说我做这件事的原因，其实我并不太关心计算模拟方面的研究，平时做软件系统比较多，一不小心接了一个水利研究单位的项目，经费其实很少，但主要想挑战一下难题（水利的研究人员很想实现），以往也没涉及过水利这个领域，所以愿意试试。后来，这个项目以开放实验室课题的方式下达到单位，名称是“水利工程力学计算遗留系统的软件再工程技术研究”，其实就是针对水利上使用的一个fortran仿真程序进行改造，使其可以用C#混合编程。

对方实验室的研究员提出一个解决方案，想让我们把这个程序用C#重写一遍。这个想法当然被我否决了，因为难以保证写出的程序和原程序输出一样，而且水利业务我不懂，没法判断是否改错了。我提出的方案是，把原有spysics程序改造成c兼容的dll链接库，把必要的接口留出，并写一个C#的适配器，这样，C#调用这个库就像使用C#对象一样简单了。我以前做过一些多语言混编的工作，可行性应该是有的。

想法是不错，关键还得动手。立项一年了，项目基本没有动，还有一年就要结项了。前两天突然接到提交中期报告的通知，有点茫然，确实不想花很大精力做这个事，项目经费实在太少，但还得做，项目不能黄掉，面子不能掉地下......

趁着这个元旦假期前的一周，没有太多事，做了做准备就开始干了。说说我的初始状态：老程序员一个，但是fortan零基础，SPHsics2D和3D仿真程序没有用过。所以第一天上手先看SPHsics2D源码带的用户文档，找了ivf和ftn95把程序跑起来，搞清楚各个.bat文件里面的细节（立项以前也看过，但不认真，没有看懂，看来世上就怕认真二字）。第二天装虚拟机环境（准备试试多种工具的配置，不想在自己电脑上搞太乱）和工具组合，失败了好几次。第三天终于配好了一个稳定且速度不错的环境，进一步研究代码结构和配置文件，晚上睡觉前提出了一个设想。第四天一早就开始迫不及待的开始整理新的结构、改一些代码，
到晚上12：30，终于长吁一口气，自己在中期检查报告里写的进度终于搞定了，而且剩下的一半工作没有难度了，只剩工作量！

于是今天写下来，以免时间太长忘得一干二净。

---------------------------------------------------------------------------------------------------------------------------------------------------------------

1、源码版本、工具选型

1）SPHYsics源码：

开源网站地址：https://wiki.manchester.ac.uk/sphysics/index.php/Downloads
我下载的是SPHYSICS_2D_v2.2.001.zip和SPHYSICS_3D_v2.2.001.zip这个源码包，为比较新的2.2.001版本（最后更新时间是2011年1月）

2）工具选型：

虚拟机使用winXp，集成编译环境使用visual studio 2010，fortran编译器选用Intel parallel studio XE 2011（ivf2011），数据文件查看器使用paraview-3.0.2。
上述工具是试过的配合较好的组合，之前也选过win10+vs2017+ivf2018的豪华组合，无奈机器跑不动只好作罢。
安装的顺序是：先装vs2010，再安装ivf2011。这个不能错，否则ivf集成不到vs中，没法用vs编译fortran程序

2、代码重构的思路与结构设计
1）代码执行效果

SPHYsics的代码结构如下

run_directory是所有例子的运行路径，进入例子文件夹（如\run_directory\Case1）直接执行相应脚本即可产生一个。
本环境配置下的执行方法是：使用ivf2011的32位命令行工具IA-32 Visual Studio 2010 mode，用cd命令通过命令行进入这个文件夹，并执行Case1_windows_ifort.bat（因为安装的是ivf），

可得到一组生成文件，其中SPHYSICSgen_2D.exe 先生成，并由它根据case1.txt中的配置，生成其他文件例如SPHYSICS_2D.mak，接下来由SPHYSICS_2D.mak生成特定的SPHYSICS_2D.exe（从源码中选取了一部分,obj链接成exe），接下来由SPHYSICS_2D.exe根据参数生成PART_000X文件，这些是用于仿真的某一帧图像的数据文件。

将post-processing \PART2VTU_windows_ifort.bat拷贝到Case1文件夹下，运行可得到PART2VTU_2D.exe，并自动将PART_000X转换为可被paraview读取的VTUinp.pvd文件。

使用paraview读取的效果如下

2）代码分析

上述过程涉及几个文件，第一个是Case1_windows_ifort.bat

@ECHO OFF
del *.exe
del *.makset UDIRX= %CD%cd ..\..\execs\move SPHYSICSgen_2D.exe ..\execs.bakcd ..\source\SPHYSICSgen2D
del *.exenmake -f SPHYSICSgen_win_ifort.mak clean
nmake -f SPHYSICSgen_win_ifort.makIF EXIST SPHYSICSgen_2D.exe (ECHO.ECHO SPHYSICSGEN compilation Done=yesECHO.copy SPHYSICSgen_2D.exe ..\..\execs\SPHYSICSgen_2D.execd %UDIRX%copy ..\..\execs\SPHYSICSgen_2D.exe SPHYSICSgen_2D.exeSPHYSICSgen_2D.exe <Case1.txt > Case1.outcopy SPHYSICS.mak ..\..\source\SPHYSICS2D\SPHYSICS.makcd ..\..\execs\del *.objmove SPHYSICS_2D.exe ..\execs.bakcd ..\source\SPHYSICS2Ddel *.exenmake -f SPHYSICS.mak cleannmake -f SPHYSICS.makIF EXIST SPHYSICS_2D.exe (ECHO.ECHO SPHYSICScompilationDone = yesECHO.copy SPHYSICS_2D.exe ..\..\execs\SPHYSICS_2D.execd %UDIRX%copy ..\..\execs\SPHYSICS_2D.exe SPHYSICS_2D.exe SPHYSICS_2D.exe) ELSE (ECHO.ECHO SPHYSICScompilation FAILEDECHO Check you have specified the correct compiler in Case fileECHO.cd %UDIRX%)
) ELSE (ECHO SPHYSICSGEN compilation FAILEDcd %UDIRX%
)

可以看到，主要是生成了SPHYSICSgen_2D.exe，并把<Case1.txt>作为输入参数，产生Case1.out。

SPHYSICSgen_2D.exe <Case1.txt > Case1.out

进一步，还产生了SPHYSICS.mak用于生成SPHYSICSgen.exe

nmake -f SPHYSICS.mak

通过执行Case1和Case2中的脚本，会发现两个例子产生的SPHYSICSgen.exe体积是不同的。根据SPHYSICS.mak可以发现，每个Case产生的SPHYSICS.mak是不同的。下面是Case1中SPHYSICS.mak的代码

OPTIONS= /NOLOGO
COPTIONS= /03OBJFILES=energy_2D.obj recover_list_2D.obj \ini_divide_2D.obj keep_list_2D.obj \SPHYSICS_2D.obj getdata_2D.obj \check_limits_2D.obj \divide_2D.obj \movingObjects_2D.obj movingGate_2D.obj \movingPaddle_2D.obj movingWedge_2D.obj \updateNormals_2D.obj vorticity_2D.obj\periodicityCorrection_2D.obj \ac_NONE_2D.obj \kernel_correction_NC_2D.obj \ac_2D.obj \poute_2D.obj \gradients_calc_basic_2D.obj \self_BC_Dalrymple_2D.obj \celij_BC_Dalrymple_2D.obj \rigid_body_motion_2D.obj \variable_time_step_2D.obj \viscosity_artificial_2D.obj \correct_2D.obj \kernel_wendland5_2D.obj \EoS_Tait_2D.obj \densityFilter_MLS_2D.obj \ac_MLS_2D.obj \LU_decomposition_2D.obj \pre_celij_MLS_2D.obj \pre_self_MLS_2D.obj \step_predictor_corrector_2D.obj 
.f.obj:ifort $(OPTIONS) $(COPTIONS) /O3 /c $<SPHYSICS_2D.exe: $(OBJFILES)xilink /OUT:$@ $(OPTIONS) $(OBJFILES)clean:del *.mod *.obj

里面列出了构成当前SPHYSICS_2D的.obj文件，该文件将这些obj链接成可执行的.exe。对比\source\SPHYSICS2D可发现，生成的obj文件其对应的.f文件是文件夹中所有.f文件的子集。分析Case2的SPHYSICS.mak可知，构成Case2中SPHYSICS_2D.exe的是所有.f文件的另一个子集。

对照分析\source\SPHYSICSgen2D中的SPHYSICSgen_2D.f及其输入文件Case1.txt可知，SPHYSICSgen_2D.f根据Case1.txt中的选项，将相应.f文件写入SPHYSICS.mak。

Case1.txt如下所示。文件分为两列，第二列是问题，第一列是问题的答案。SPHYSICSgen_2D.f就是根据第一列进行.f文件选择的。

0  			Choose Starting options:  0=new, 1=restart, 2=new with CheckPointg, 3=restart with CheckPointing
5  			Kernel: 1=gaussian, 2=quadratic; 3=cubic; 5=Wendland
1  			Time-stepping algorithm: 1=predictor-corrector, 2=verlet, 3=symplectic, 4=Beeman
2  			Density Filter: 0=none, 1=Shepard filter, 2=MLS
30                      ndt_FilterPerform ?
0  			Kernel correction 0=None, 1=Kernel correction, 2=Gradient kernel Correction
1   			Viscosity treatment: 1=artificial; 2=laminar; 3=laminar + SPS
0.3   			Viscosity value( if visc.treatment=1 it's alpha, if not kinem. visc approx 1.e-6)
0                       vorticity printing ? (1=yes)
1    			Equation of State: 1=Tait's equation, 2=Ideal Gas, 3= Morris
2     			Maximum Depth (h_SWL) to calculate B
10     			Coefficient of speed of sound (recommended 10 - 40 ) ??
2    			Boundary Conditions: 1=Repulsive Force; 2=Dalrymple
15                      ndt_DBCPerform ? (1 means no correction)
1  			Geometry of the zone:  1=BOX, 2=BEACH, 3=COMPLEX GEOMETRY
2  			Initial Fluid Particle Structure: 1= SC, 2= BCC
4.0,4.0      		Box dimension LX,LZ?
0.03,0.03   	    	Spacing dx,dz?
0    			Inclination of floor in X ( beta ) ??
0,0,0                   Periodic Lateral boundaries in X, Y, & Z-Directions ? (1=yes)
0    			Add wall
0     			Add obstacle (1=y)
0     			Add wavemaker (1=y)
0     			Add gate (1=y)
0     			Add Floating Body (1=yes)
2      			Initial conditions: 2) particles on a staggered grid 
0      			Correct pressure at boundaries ?? (1=y)
0.03,1.     		Cube containing particles :  XMin, Xmax ??
0.03,2.    	        Cube containing particles :  ZMin, Zmax ??
0     			Fill a new region
3.0,0.02      	        Input the tmax and out
0.   			initial time of outputting general data
0.0005,1.0,-1.0    	For detailed recording during RUN: out_detail, start, end
0.0001,1     		Input dt?? , i_var_dt ??
0.2    			CFL number (0.1-0.5)
0.92     		h=coefficient*sqrt(dx*dx+dz*dz):  coefficient ???
0   			Use of Riemann Solver: 0=None, 1=Conservative (Vila), 2=NonConservative (Parshikov)
3  			Which compiler is desired: 1=gfortran, 2=ifort, 3=win_ifort, 4=Silverfrost FTN95
1  			Precision of XYZ Variables: 1=Single, 2=Double

可以对照一下SPHYSICSgen_2D.f的一段代码，其中"read(*,*) i_restartRun"就是针对Case1.txt中的第一行数据进行读取。

      write(*,*) 'Choose Starting options:   Start new RUN = 0 ' write(*,*) '                       : Restart old RUN = 1 'write(*,*) '     with CheckPointing:   Start new RUN = 2 'write(*,*) '                       : Restart old RUN = 3 'read(*,*) i_restartRunwrite(*,*) i_restartRunC     KERNELwrite(*,*)'Choose KERNEL'write(*,*)'Gaussian        = 1'write(*,*)'Quadratic       = 2'write(*,*)'Cubic- Spline   = 3'write(*,*)'Quintic Wendland= 5'read(*,*) i_kernelwrite(*,*) i_kernel

由此，大致可得出结论：

首先生成SPHYSICSgen_2D.exe，再根据Case1.txt生成SPHYSICS_2D.exe，然后运行SPHYSICS_2D.exe，根据INDAT文件生成仿真数据PART_000X。

3）代码重构思路

由于目前SPHYSICSgen_2D.exe和SPHYSICS_2D.exe是两个独立的程序，且SPHYSICS_2D还是一个由Case.txt配置生成的，执行特定任务的功能子集，因此考虑将SPHYSICSgen_2D和SPHYSICS_2D的代码打包成一个动态链接库dll，根据需要给出相应函数入口，由C#等程序调用。

后续使用这个库将会比较方便，不需要每次生成新的SPHYSICS_2D，需要哪些函数都可以从dll中调取，而执行特定任务可以通过写Case.txt文件来做到命令的批处理。

但是，进一步分析\source\SPHYSICS2D中的文件会发现，存在多个不同的文件函数名称相同的情况，直接引入工程中编译，会报同名函数的错误。

分析SPHYSICSgen_2D.f代码会发现，有多个类似下图的代码。这些例如i_kernelcorrection的变量都是从前面"read(*,*) i_kernelcorrection" 中的语句获得值。也就是说，Case1.txt指定了一些配置，i_kernelcorrection的值不同，对应到生成SPHYSICS_2D时会选择不同的.f文件引入。打开这些文件，如“ac_kgc_2d.f”、“ac_kc_2d.f”会发现其函数名都是ac_main。

      !- Kernel Correctionsif (i_kernelcorrection.eq.0) thenwrite(22,FMT1)TAB,'ac_NONE_2D.o \'write(22,FMT1)TAB,'kernel_correction_NC_2D.o \'elseif (i_kernelcorrection.eq.2) thenwrite(22,FMT1)TAB,'ac_KGC_2D.o \'write(22,FMT1)TAB,'kernel_correction_KGC_2D.o \'write(22,FMT1)TAB,'pre_self_KGC_2D.o \'write(22,FMT1)TAB,'pre_celij_KGC_2D.o \'elseif (i_kernelcorrection.eq.1) thenwrite(22,FMT1)TAB,'ac_KC_2D.o \'write(22,FMT1)TAB,'kernel_correction_KC_2D.o \'endif

于是我猜测，应该是在生成SPHYSICS_2D时会引入多个子程序或函数，其中一些子程序或函数有多个可供的替代函数，在每一个生成的SPHYSICS_2D里，不会出现这些替代函数的冲突（若冲突，则不会链接出exe）。

这让我想到“多态”的情况，就是同一个操作泛型，不同的实现方法。后续对代码进行排查进一步证实了这个想法。我排查的方法很简单，就是让编译器先编译报错，我在debug的时候注意观察这些同名的方法。不过，fortran的编译查错真的很痛苦，它通常不会报出行号。

于是我准备试一试，通过手工构造一些函数达到多态效果。其实我对fortran没有基础，也没有仔细看过教程，好在是个老程序员，一边靠着经验，一边靠着百度，就开始改造了。

SPHYSICS的代码是fortran77格式的，比较难调。据说fortran可以面向对象，但这个77版大概是面向不了对象啦，我就把它当c程序或者matlab的过程来改吧。

我在vs里选择ivf的fortran 模板，新建了一个Dynamic library 项目，把SPHYSICS2D的代码都贴了进去。

拿到的SPHYSICS是所有文件都在一个文件夹里面，为了便于修改和对照，构造了以下结构，将函数名相同的文件放到一个文件夹下，文件夹名为他们的同名函数名。

例如上图中的ac_2D.f 和 ac_Conservative_2D.f，他们的部分原始代码分别如下，有一部分相同，但又有一部分不同。

     subroutine ac_main
cinclude 'common.2D'c
c  ...  store useful arrays
c!- Need to zero each object for multiobjects -bigUdot   = 0.0bigWdot   = 0.0X_Friction   = 0.0Z_Friction   = 0.0nb_inFriction = 0

      subroutine ac_main
cinclude 'common.2D'c
c  ...  store useful arrays
c!- Need to zero each object for multiobjects -bigUdot   = 0.0bigWdot   = 0.0X_Friction   = 0.0Z_Friction   = 0.0nb_inFriction = 0do i=nbfm+1,nbc        -- Zeroing Variables for Free-Moving Objects --

在ac_main文件夹中建立ac_main.f文件（在下一步加入配置项时由该文件进行路由，选择合适文件），作为该文件夹中所有函数文件的泛型，并将这些文件的函数名改为与其文件名一致的名称，如ac_2D.f中原有的ac_main函数改为ac_2D，ac_Conservative_2D.f中原有ac_main函数改为ac_Conservative_2D。

对所有文件夹都做如此的操作，编译，并解决一些小问题（如参数列表不同，则无需做泛型文件，直接在合适的位置修改）。

以这种方法进行修改和编译，最终可以生成dll。

3、重要代码段修改说明

但这还达不到效果，因为原有SPHYSICSgen_2D.exe和SPHYSICS_2D.exe配合生成数据的方式中，是通过配置文件Case1.txt进行函数选择的，把exe改造成dll后失去了自动读取文件的能力。

因此，在dll的工程文件中加入SPHYSICSgen_2D.f文件，并改造之：

cc---------sw--------------------------------------------------------------      open (2, file='case.txt', status='old')
cc---------sw--------------------------------------------------------------       write(*,*) 'Choose Starting options:   Start new RUN = 0 ' write(*,*) '                       : Restart old RUN = 1 'write(*,*) '     with CheckPointing:   Start new RUN = 2 'write(*,*) '                       : Restart old RUN = 3 '
cc---------sw-------------------------------------------------------------- read(2,*) i_restartRun
cc---------sw--------------------------------------------------------------

加入 open (2, file='case.txt', status='old')，读取配置文件，并将原有read(*,*) i_restartRun标准输入改为read(2,*) i_restartRun，读取2号文件中的值。相应的修改很多，都是这个原理。

此外，由于不需要生成SPHYSICS_2D文件了，而应该生成构成SPHYSICS_2D文件的函数列表，因此找到subroutine tocompile_win_ifort子程序，将其中的

      open(22,file='SPHYSICS.mak')write(22,FMT) 'OPTIONS= /NOLOGO'write(22,FMT) 'COPTIONS= /03'write(22,FMT)write(22,FMT) 'OBJFILES=energy_2D.obj recover_list_2D.obj \'write(22,FMT1)TAB,'ini_divide_2D.obj keep_list_2D.obj \'      write(22,FMT1)TAB,'SPHYSICS_2D.obj getdata_2D.obj \'write(22,FMT1)TAB,'check_limits_2D.obj \'write(22,FMT1)TAB,'divide_2D.obj \'write(22,FMT1)TAB,'movingObjects_2D.obj movingGate_2D.obj \'write(22,FMT1)TAB,'movingPaddle_2D.obj movingWedge_2D.obj \'write(22,FMT1)TAB,'updateNormals_2D.obj vorticity_2D.obj\'write(22,FMT1)TAB,'periodicityCorrection_2D.obj \'

改为如下

      open(22,file='SPHYSICS.fun')write(22,FMT) 'energy_2D'write(22,FMT) 'recover_list_2D'write(22,FMT) 'ini_divide_2D'write(22,FMT) 'keep_list_2D'write(22,FMT) 'SPHYSICS_2D'write(22,FMT) 'getdata_2D'write(22,FMT) 'check_limits_2D'write(22,FMT) 'divide_2D'write(22,FMT) 'movingObjects_2D'write(22,FMT) 'movingGate_2D'write(22,FMT) 'movingPaddle_2D'write(22,FMT) 'movingWedge_2D'write(22,FMT) 'updateNormals_2D'write(22,FMT) 'vorticity_2D'write(22,FMT) 'periodicityCorrection_2D'

目的是输出一个.fun文件，便于读取应该配置哪些同名函数实现“多态”效果。

接下来定义一个sph_module.f文件，在里面读取.fun文件并为一组全局变量赋值，这一组全局变量就是后续要实现“多态”，进行文件选择的关键。从下图可以看出，.fun文件被读出，并赋值给相应全局变量。

        open(3,file='SPHYSICS.fun',status='old')do 10 i=1,32000read(3,*,iostat=stat1) contentif (stat1.ne.0) goto 10c-------------ac_main_str-------------------------------------------------------            trimStr=trim(content)if(trimStr.eq."ac_2D") thenac_main_str=trimStrwrite(*,*) ac_main_strendifif(trimStr.eq."ac_Conservative_2D") thenac_main_str=trimStrwrite(*,*) ac_main_strendifif(trimStr.eq."ac_KC_2D") thenac_main_str=trimStrwrite(*,*) ac_main_strendifif(trimStr.eq."ac_KGC_2D") thenac_main_str=trimStrwrite(*,*) ac_main_strendif
c-------------celij_str-------------------------------------------------------

在上一节中建立的ac_main.f等“多态”路由文件还没有真正编写内容，现在就要起作用啦！

对ac_main.f等“多态”路由文件编写如下代码，可以看出，当调用ac_main函数时，会根据配置文件赋值的全局变量调用相应“替代”函数，从而达到“多态”的效果。

      subroutine ac_mainuse sph_moduleinclude 'common.2D'if(ac_main_str.eq."ac_2D") thencall ac_2Dendifif(ac_main_str.eq."ac_Conservative_2D") thencall ac_Conservative_2Dendifif(ac_main_str.eq."ac_KC_2D") thencall ac_Conservative_2Dendifif(ac_main_str.eq."ac_KGC_2D") thencall ac_KGC_2Dendif  returnend

在SPHYSICSgen_2D.exe和SPHYSICS_2D.exe配合生成数据的方式中，是通过选择.obj文件，生成无重名函数的SPHYSICS_2D.exe。而通过上述方法，可以在dll中根据配置文件调用相应方法，从而实现无需每次编译，也达到dll可重用的效果。

4、动态库的生成与C#混编效果

用上述方法修改一部分代码的函数头部，将其声明为DLLEXPORT，便于以编译C的方式访问，编译通过后生成dll。

      subroutine SPHYSICSgen_2D!DEC$ ATTRIBUTES DLLEXPORT::SPHYSICSgen_2D!DEC$ ATTRIBUTES STDCALL,ALIAS:'SPHYSICSgen_2D'::SPHYSICSgen_2D

将dll加入一个C#的控制台工程，为dll写一个适配类FortranMethod，就可以像C#类那样使用SPHYSICS了。

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Runtime.InteropServices; namespace testSPH
{class Program{static void Main(string[] args){FortranMethod.SPHYSICSgen_2D();//FortranMethod.SPHysics();}}public static class FortranMethod{[DllImport("SPH2D.dll")]public static extern void SPHYSICSgen_2D();//[DllImport("SPH2D.dll")]//public static extern void SPHysics();}
}

运行的结果和直接运行bat文件一样。不过在C#下调用，会比用fortran的原生exe慢不少。