写在前面——如何阅读OpenPhase源码
OpenPhase的源码真的是挺难读的，它是一整套完全的很成熟的求解相场模型的框架，其模块化完成度很高(尤其设计的数据结构面面俱到，号称能处理所有类型的数据)，所以用它来解方程时写的代码很少，但带来一个问题是你必须先读懂它：高度集成的代码让初学者“一进去就出不来”(从而迷失在各个类中)或“压根进不去”(一行代码就完成一个大功能，根本不知道怎么入手)。
我的建议是：
- 善用Doxygen
OpenPhase支持Doxygen，其注释也较全面，使用Doxygen生成类之间的关系图，方便理清思路。
附录有Doxygen生成的dot图的涵义说明。 - 使用gdb调试
刚开始可能不知道OpenPhase里各个参数的值是多少，可以使用gdb逐步调试且输出关键参数的值，增加理性认识。也可以在源码中加cout来输出，但这样有些麻烦。
本文尝试对HeatEquation这个算例进行详细解析，试图理清OpenPhase的计算思路。
模拟参数
本例模拟参数分布在四个opi文件中：
ProjectInput.opi文件
此文件中的参数是通过Settings类的ReadInput函数读入。
针对本例，读入的重要参数有：
- 晶粒个数是2: Ngrans=2
- 相的个数是2: Nphses=2
- 格点数目为：Nx=64, Ny=1, Nz=64
- 网格间距为：dx=1.0
- 界面宽度占格点的数目为：iWidth=0
- 系统实际温度：T=500
- OpenMP线程个数：nOMP=2
BoundaryConditions.opi文件
此文件中的参数通过BoundaryConditions类的ReadInput函数读入。
- 六个边界的边界条件：BC0X=Free, BCNX=Free, BC0Y=Free, BCNY=Free, BC0Z=Free, BCNZ=Free
Temperature.opi文件
此文件中的参数通过Temperature类的ReadInput函数读入。
- 初始温度梯度所在位置的坐标：R0X=32, R0Y=32, R0Z=32
- 初始温度梯度：DT_DRX=0, DT_DRY=0, DT_DRZ=0
- 冷却速度：DT_Dt=0
Heat.opi文件
此文件中的参数通过Heat类的ReadInput函数读入。
- 相0的热扩散率：ThermalDiffusivity_0=13.0e-5
- 相0的热容：HeatCapacity_0=0.465e-3
- 相1的热扩散率：ThermalDiffusivity_0=0.0
- 相0的热容：HeatCapacity_0=5.0
初始化
相场的初始化：
- 依次读入三个维度的格点数目、网格间距、界面宽度、相的个数
- 分配Fields场空间
Fields是Storage3D
Size_X=Nx+2=66, Size_Y=Ny+2=3, Size_Z=Nz+2=66
同时开辟出一块大小为Size_XSize_YSize_Z=13068的Node类型的内存，将初始地址赋给指针Array。 - 分配Fractions场空间
Fractions是Storage3D
Size_X=Nx+2=66, Size_Y=Ny+2=3, Size_Z=Nz+2=66, Size_N=Nphases=2
但它开辟出一块的double类型的内存大小为Size_X*Size_Y*Size_Z*Size_N=26136，同时将初始地址赋给指针Array。 - 分配Laplacians场空间
Laplacians场的类型及操作与Fields场相同。 - 分配Gradients场空间
Gradients场的类型是Storage3D - 分配Normals场空间
Normals同Gradients场。 - 分配Flags场空间
Flags场是Storage3D - 分配FieldIndex场空间
FieldIndex同Flags场。 - 分配FieldsStatistics场空间
FieldsStatistics场是GrainInfo类型的场，其Allocate函数传入晶粒个数Ngrains和相个数Nphases，此处初始化就直接将Ngrains设为1，该场的成员变量GrainStorage是一个元素为Grain类型的容器，Allocate函数将该容器的大小(即里面元素的个数)设定为Ngrains*Nphases=1*2=2，同时对该场的成员变量Nphases赋值。 - 对相场的成员变量TotFractions(是元素为double的容器)的尺寸设定为Nphases=2，传入iWidth，对RefVolume赋值，根据Nx等与iWidth的关系判断，此处为0。
温度场的初始化
与相场的初始化的逻辑相同，也是读入网格格点、相的个数、网格间距。对场Tx、Txdx、qdot三个场分配空间，注意Txdx的Array指向dVector3类型的变量。初始温度T0=500, 初始温度梯度dT_dr是一个dVector3的变量，此处调用它的set_to_zero成员函数将其成员变量Storage[]清零(这里调用了c++的memset函数)。初始温度梯度所在的位置r0也是一个dVector3，其set_to_zero同理。
热传导方程求解器的初始化
与以上的逻辑相同，相继读入网格格点、相的个数、网格间距。成员变量PhaseThermalDiffusivity是一个类型为Storage
边界条件的初始化
其Initialize函数传入的是Settings类，但此处初始化中直接将其六个int型的成员变量BC0X、BCNX、BC0Y、BCNY、BC0Z、BCNZ都用OpenPhase的整型常量Periodic=0来赋值。
以上可以看出，OpenPhase在设计Initialize函数时并没有完全真的读入相应的参数并初始化，而是采取类似于定义一个变量并赋无意义的初值这样的做法。下面是真正的读取输入文件。
边界条件读取输入文件
调用ReadInput函数，传入的形参是FileName字符串，实参是ProjectDir+BCInputFileName，这两个字符串是命令空间OpenPhase的全局变量，默认是ProjectDir = “ProjectInput/“和BCInputFileName = “BoundaryConditions.opi”。
调用Type.h定义的ReadParameterF函数读取.opi中的数据(详见后面的解析)，取得了该文件的字符串，然后通过判断字符串，将相应的边界条件的整型数值传入BCOX、BCNX等六个成员变量：0是周期性边界；1是非流动边界；2是自由边界；3是固定边界。
温度场读取输入文件
通过ReadInput函数读取初始温度梯度所在位置的坐标，并赋值给其成员变量r0，即r0[0]=R0X=32,r0[1]=R0Y=32,r0[2]=R0Z=32，起始r0[]是一个运算符重载，实际将参数传给r0的成员变量Storage数组中。同理将初始温度梯度赋值给dT_dr，即dT_dr[0]=DT_DRX=0,dT_dr[1]=DT_DRY=0,dT_dr[2]=DT_DRZ=0。最后将冷却速度赋值给成员变量dT_dt。
热传导求解器读取输入文件
Heat的ReadInput函数有些特殊，其会对不同的相进行循环，对于相的指标n，分别将热导率和热容存入PhaseThermalDiffusivity和PhaseHeatCapacity中，注意[]依然是一个运算符重载，实际是将数值存入变量的Array指针指向的地址中。
经过上面的初始化，各个物理场中存储了正确的模拟参数，下面是初始条件设置。
相的分布设定
Initializations::Single函数
该函数首先调用PhaseField类的FieldStatistics成员变量的add_grain函数，传入的参数是PhaseIndex这个重要参数，此处的实参是0，即代表0号相。注意GrainStorage在之前相场的初始化分配FieldsStatistics时已经设定了，其size()=2，经过add_grain后，返回locIndex=0，GrainStorage[0].Exist=true。
然后对格点进行循环，对Flags、Fields、FieldIndex这三个成员变量进行设定。注意到之前初始化时Allocate后开辟的各自的Array指针的地址都是Size_X*Size_Y*Size_Z=66*3*66=13068大小，故循环是for(int i = 0; i < Nx+2; ++i)这样的类型。注意赋值时使用Storage
(Size_Y*x+y)*Size_Z+z=(Size_Y*(Size_X-1)+Size_Y-1)*Size_Z+Size_Z-1=Size_X*Size_Y*Size_Z-1
就这样取得Array指针指向的所有地址后，再对其赋值。对于Flags(i,j,k)，其返回值是int型，直接赋值0。而对于Fields，其返回值是Node类型，其再调用Node类型的set函数(注意是有两个形参的set)。具体的set过程为：Node类的Fields成员变量(注意与前面的Fields区别)是一个元素为FieldEntry的容器，FieldEntry是一个结构体，从Fields.begin()开始做循环时，发现里面根本没有东西，就新建一个NewEntry，并将传入的两个参数,一个是locIndex=0赋给NewEntry的index，一个是1.0赋给NewEntry的value，NewEntry的另外两个值为0，然后将NewEntry压入Fields中。后面对于FieldIndex(i,j,k)就很自然地取得刚刚的index。
最后的Finalize函数由多个部分组成：
- FinalizeStepOne(BoundaryCondition&)函数。
它使用了OpenMP的并行计算方法，这个地方不好用gdb调试，所以可以直接在源码中用cout输出。首先用omp_get_num_threads()得到线程总数，此处为nOMP=2，然后用omp_get_thread_num()得到线程标识，分别为0和1，然后将x方向的格点分为两段，分别是从1到32和从33到64，然后对三维的格点进行循环，判断Flags(i,j,k)是否为0，因为上面已经设定此值为0，所以此处不进入if处理。 - SetBoundaryConditions(BoundaryCondition&)函数。
此函数传入BoundaryCondition类的对象，实际是接着调用BC的SetX、SetY、SetZ三个成员函数，注意这里有个函数重载，传入的参数分别是PhaseField类的Fields和FieldIndex成员变量。SetX函数是根据不同的边界条件对Fields和FieldIndex进行设定。如对Periodic边界条件，让0号和64条件相同，1和65号相同。注意这里面设定的对应编号。对于左边界，如果是NoFlux边界条件，则让0号和1号相同，如果是Free边界，则让0号+2号等于2倍的1号。Fixed边界则不做处理。对于右边界，有同样的处理。SetY和SetZ亦同。 - SetFlags函数
其会对Interface(i,j,k)函数作判断，然后满足某些条件后将Flags(i,j,k)=Interface(i,j,k)+1。此处没有满足if条件，所以不做处理。 - SetFlagsBC函数
此处调用SetX等对Flags进行边界条件的处理。 - CalculateFractions函数
计算相分数，此处也会作一些判断，但此时仍然未满足条件，只是将所有地方的Fractions都置为1。 - CalculateLaplacians函数
计算拉普拉斯算子，仍然没满足条件，所以未分析。 - CalculateGradients函数
计算梯度，处理同上。 - SetGradientsBC函数
设置梯度边界条件 - CalculateNormals函数
计算法线向量 - CalculateVolumes函数
计算相体积。
Initializations:Sphere函数
Sphere函数与Single函数逻辑类似，Single相当于Sphere的一个初级版本，很多变量都是0，因此对相场的设置很简单，而Sphere则进一步将不同的相分开，赋予不同的属性，因此更能加深对代码的理解。
这部分的具体步骤不再详述。
对条件的进一步设定
此部分对于温度场有了进一步设置，如初始温度Tx.T0和热源Tx.qdot。
求解温度场
主要利用Heat类的Solve成员函数。
判断稳定性
判断时间步长与网格间距之间的关系是否满足稳定性要求。
这里设定为dt/(dx*dx)<0.5。
Temperature.SetBoundaryConditions(BC)函数
传入Tx参数对Ghost Nodes设定边界条件。
SetEffectiveProperties(Phase)函数
这一步是将相场的热容、热扩散率、密度等材料物性参数传递到Heat类的对应的成员变量中。
CalculateLaplacian(Temperature)函数
这一步使用有限差分方法计算拉普拉斯项——二阶微分项，具体的差分格式是中心差分。
UpdateTemperature(Temperature, dt)函数
求解热传导方程,更新温度场
输出
相场的输出用Phi.WriteVTK(tStep)函数，格式是标准的VTK格式，格式的具体解析见此文。
具体输出的内容为：
(1)格点的编号，共Nx*Ny*Nz个点，排列次序为：
1 | 1 | 1
2 | 1 | 1
…| … | …
64 | 1 | 64
(2)输出第一个场——Interfaces，1代表界面，0代表体相
(3)输出第二个场——PhaseFields，
(4)输出第三个场——PhaseFraction，该场的个数与相的个数Nphases相等
(5)输出第四个场——ColorScale，
(6)输出第五个场——Flags，
注意事项
- 不同的变量用的编号不同，有的是从0号到63号，有的是从1号到64号，还有65号的，这期间应涉及一个特殊处理——Ghost nodes，暂时没弄明白怎么回事，留待以后解决。
附录一：Doxygen的dot图的涵义
上图中的box的涵义：
- 实心的灰色box代表绘制此图所基于的结构体或类
- 黑色边界的box意味着此结构体或类有记录
- 灰色边界的box意味着此结构体或类无记录
- 红色边界的box意味着此有记录的结构体或类的关系没有被完全绘制出来
上图中的箭头的涵义：
- 深蓝色箭头表示两类之间的public继承关系(上方的类是父类，下方的为子类)
- 深绿色箭头表示两类之间的protected继承关系
- 深红色箭头表示private继承关系
- 紫色虚线箭头表示其所指的类或结构体被包含或使用，上面的label是使用此类的变量
- 黄色虚线箭头表示模板实例和模板类之间的关系，上面的label是该实例的模板参数
附录二：ReadParameterF函数
ReadParameterF函数有四个形参：文件流inp、字符串Key、整型变量necessary、字符串defaultVal。传入的实参是BoundaryCondition.opi，将其赋值给inp，将字符串BCOX、BCNX等赋值给Key。
开始读入文件时，截取文件开头一直到字符的内容，赋值给tmp，如果没有这个字符，则抛出一个错误信息；然后再继续读取文件，将空格之前的字符串赋值给ReadKey字符串，如BC0X，将ReadKey与Key通过compare函数对比，如果两者不同的话，则继续读取；如果相同的话，再将后续到字符:的内容赋值给tmp，接着将直到回车符的字符串赋值给tFileName，如Free。为了使格式整齐，使用字符串的erase等函数删除空格，然后调用Info类的GetStandard输出。将tFileName填入ReturnValue中，返回此值。
