使用PRISMS-PF模拟二维和三维枝晶生长

引子

本例将使用PRISMS-PF模拟凝固过程中的枝晶生长。因为原先的PRISMS-PF架构中在构建右端项时只有接收一步解,而枝晶的浓度场控制方程与相场参量的变化率相关,即需要用到相场的前后两步解,所以需要扩展PRISMS-PF的原先架构,本例也主要是说明怎样对原程序进行扩展。

获得源码

源码在这里

扩展过程

定义要使用解的个数

在parameter.h中新增宏定义numSolution,其决定方程中用到几步解,比如在枝晶生长的模拟中用到两步解,其就等于2,之前的所有模拟中都等于1。

1
#define numSolution 2

这里涉及头文件的包含顺序,放在后面定义的话,会不认识。

新增前一步解和临界解的存储对象

仿照solutionSet的定义,在相同位置(具体就是MatrixFreePDE类声明)增加oldSolutionSet和tempSolutionSet,同时初始化空间,并且根据初始条件初始化值。
注意这个地方初始化时不要跟当前解使用同一个地址,否则它们自动跟着当前解变化,应该另外开辟内存地址:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
vectorType *U,*oldU, *tempU, *R;
if (iter==0){
U=new vectorType; R=new vectorType;
oldU=new vectorType; tempU = new vectorType;
solutionSet.push_back(U); residualSet.push_back(R);
oldSolutionSet.push_back(oldU);
tempSolutionSet.push_back(tempU);
matrixFreeObject.initialize_dof_vector(*R, fieldIndex); *R=0;
}
else{
U=solutionSet.at(fieldIndex);
oldU=oldSolutionSet.at(fieldIndex);
tempU=tempSolutionSet.at(fieldIndex);
}
matrixFreeObject.initialize_dof_vector(*U, fieldIndex); *U=0;
matrixFreeObject.initialize_dof_vector(*oldU, fieldIndex); *oldU=0;
matrixFreeObject.initialize_dof_vector(*tempU, fieldIndex); *tempU=0;

更新解

每次计算开始前,将solutionSet暂存在tempSolutionSet中,计算完毕后,使用tempSolutionSet对oldSolutionSet赋值(可以在基类中增加这么一个更新解的成员函数)。
注意,用当前步的解给上步解赋值时,不能直接使用等号,像这样:

1
oldSolutionSet = solutionSet;

这样会将当前解的地址直接赋给上步解,导致两个解完全同步。
需要只是传递值:

1
2
3
4
5
for(unsigned int fieldIndex=0; fieldIndex<fields.size(); fieldIndex++){
for (unsigned int dof=0; dof<solutionSet[fieldIndex]->local_size(); ++dof){
oldSolutionSet[fieldIndex]->local_element(dof) = solutionSet[fieldIndex]->local_element(dof);
}
}

将计算右端项的输入改为盛放两步解的容器并赋值

(1)将computeRHS和getRHS以及residualRHS的输入参数src改为容器,注意需要在多个地方修改:matrixFreePDE原型修改、coupled的函数定义实现、coupled_function中函数实现修改
(2)需要在计算右端项之前,对src容器赋值:通过numSolution判断读入多少个解
这个将原来盛放一步解的变量改成盛放此变量的容器,涉及多个函数及多个地方的修改,一定要仔细。
在getRHS中添加能读取上一步解的功能时,注意学习step48。两者对应关系就是:

1
2
scalar_vars VS current
old_scalar_vars VS old

在residualRHS中将输入变成容器时,注意将容器中的元素设置为存储指针,这样就直接指向原来的两个解:

1
2
3
4
std::vector<std::vector<modelVariable<dim>>*> modelVarListList;
modelVarListList.push_back(&modelVarList);
if(numSolution == 2)
modelVarListList.push_back(&oldModelVarList);

计算结果

二维结果:

三维结果

后记

(1)书写方程各项时一定注意,将其弱形式仔细写下来,尤其是分部积分时对拉普拉斯算子分解时的符号变化。
(2)学会使用天山折梅手,化简为繁;学会使用奥卡姆剃刀,化繁为简。
(3)时间步长和网格间距会导致收敛性
(4)求解精度也会影响枝晶形貌
(5)有限单元的插值方式会直接决定hessian矩阵的计算