想写的原因很简单，怕时间长了，自己也不记得过程和方法了。正好可以写下来，以后和有这方面需要的朋友们探讨。
首先说说我做这件事的原因，其实我并不太关心计算模拟方面的研究，平时做软件系统比较多，一不小心接了一个水利研究单位的项目，经费其实很少，但主要想挑战一下难题（水利的研究人员很想实现），以往也没涉及过水利这个领域，所以愿意试试。后来，这个项目以开放实验室课题的方式下达到单位，名称是“水利工程力学计算遗留系统的软件再工程技术研究”，其实就是针对水利上使用的一个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 *.mak
set UDIRX= %CD%
cd ..\..\execs\
move SPHYSICSgen_2D.exe ..\execs.bak
cd ..\source\SPHYSICSgen2D
del *.exe
nmake -f SPHYSICSgen_win_ifort.mak clean
nmake -f SPHYSICSgen_win_ifort.mak
IF EXIST SPHYSICSgen_2D.exe (
ECHO.
ECHO SPHYSICSGEN compilation Done=yes
ECHO.
copy SPHYSICSgen_2D.exe ..\..\execs\SPHYSICSgen_2D.exe
cd %UDIRX%
copy ..\..\execs\SPHYSICSgen_2D.exe SPHYSICSgen_2D.exe
SPHYSICSgen_2D.exe <Case1.txt > Case1.out
copy SPHYSICS.mak ..\..\source\SPHYSICS2D\SPHYSICS.mak
cd ..\..\execs\
del *.obj
move SPHYSICS_2D.exe ..\execs.bak
cd ..\source\SPHYSICS2D
del *.exe
nmake -f SPHYSICS.mak clean
nmake -f SPHYSICS.mak
IF EXIST SPHYSICS_2D.exe (
ECHO.
ECHO SPHYSICScompilationDone = yes
ECHO.
copy SPHYSICS_2D.exe ..\..\execs\SPHYSICS_2D.exe
cd %UDIRX%
copy ..\..\execs\SPHYSICS_2D.exe SPHYSICS_2D.exe
SPHYSICS_2D.exe
) ELSE (
ECHO.
ECHO SPHYSICScompilation FAILED
ECHO Check you have specified the correct compiler in Case file
ECHO.
cd %UDIRX%
)
) ELSE (
ECHO SPHYSICSGEN compilation FAILED
cd %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= /03
OBJFILES=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_restartRun
write(*,*) i_restartRun
C KERNEL
write(*,*)'Choose KERNEL'
write(*,*)'Gaussian = 1'
write(*,*)'Quadratic = 2'
write(*,*)'Cubic- Spline = 3'
write(*,*)'Quintic Wendland= 5'
read(*,*) i_kernel
write(*,*) 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 Corrections
if (i_kernelcorrection.eq.0) then
write(22,FMT1)TAB,'ac_NONE_2D.o \'
write(22,FMT1)TAB,'kernel_correction_NC_2D.o \'
elseif (i_kernelcorrection.eq.2) then
write(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) then
write(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
c
include 'common.2D'
c
c ... store useful arrays
c
!- Need to zero each object for multiobjects -
bigUdot = 0.0
bigWdot = 0.0
X_Friction = 0.0
Z_Friction = 0.0
nb_inFriction = 0subroutine ac_main
c
include 'common.2D'
c
c ... store useful arrays
c
!- Need to zero each object for multiobjects -
bigUdot = 0.0
bigWdot = 0.0
X_Friction = 0.0
Z_Friction = 0.0
nb_inFriction = 0
do i=nbfm+1,nb
c -- 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,32000
read(3,*,iostat=stat1) content
if (stat1.ne.0) goto 10
c-------------ac_main_str-------------------------------------------------------
trimStr=trim(content)
if(trimStr.eq."ac_2D") then
ac_main_str=trimStr
write(*,*) ac_main_str
endif
if(trimStr.eq."ac_Conservative_2D") then
ac_main_str=trimStr
write(*,*) ac_main_str
endif
if(trimStr.eq."ac_KC_2D") then
ac_main_str=trimStr
write(*,*) ac_main_str
endif
if(trimStr.eq."ac_KGC_2D") then
ac_main_str=trimStr
write(*,*) ac_main_str
endif
c-------------celij_str-------------------------------------------------------在上一节中建立的ac_main.f等“多态”路由文件还没有真正编写内容，现在就要起作用啦！
对ac_main.f等“多态”路由文件编写如下代码，可以看出，当调用ac_main函数时，会根据配置文件赋值的全局变量调用相应“替代”函数，从而达到“多态”的效果。
subroutine ac_main
use sph_module
include 'common.2D'
if(ac_main_str.eq."ac_2D") then
call ac_2D
endif
if(ac_main_str.eq."ac_Conservative_2D") then
call ac_Conservative_2D
endif
if(ac_main_str.eq."ac_KC_2D") then
call ac_Conservative_2D
endif
if(ac_main_str.eq."ac_KGC_2D") then
call ac_KGC_2D
endif
return
end在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慢不少。
更进一步，如果想使用SPHYSICS的其他函数，可以在这些函数上加上DLL导出标志。
5、参考资料
作为fortran小白，能把这个改造做完实在是不容易，参考了很多牛人的方法，在此表示感谢。
仅列出最重要的博客：
http://fcode.cn/guide-30-1.html
