MOOSE相场模块的内核模型

引子

moose的相场模块提供了通用的求解相场模型的算法,其通式采用的是自由能泛函的原始形式,只要用户要求解的模型满足这里内置的方程形式,那么用户仅需要提供自由能的导数和材料参数,就可以迅速进行模拟。比如教程中的调幅分解模型。
如果用户不能知道相场的原始形式,也可以自己开发模型,比如用于模拟枝晶生长的Kobayashi模型,就不满足上述规则,此时可以自己创建模型。

牛顿迭代法

moose内核模型的表达形式是将原来的控制方程中的右端项都移动到左端,得到这样的形式:
$$
R_i(u)=0, \qquad i=1,\ldots, N
$$
其中,N是分量的个数。
如果写成有限元函数的形式:
$$
R_i(u_h)=0, \qquad i=1,\ldots, N
$$
其中:
$$
u \approx u_h = \sum_{j=1}^N u_j \phi_j
$$
即,原来的连续的场变量由节点上的离散的系数来代替。

单变量求根

现在仅考虑单变量情形:
$$
f(x)=0
$$
moose采用Newton迭代法来对其求根,该方法的基本思想是:将非线性方程的求根问题,转化成某个线性方程的求根。比如,对该方程的左端项进行泰勒展开,且只取其线性主部:
$$
0=f(x)=f(x_0)+f’(x_0)(x-x_0)
$$
那么,解就是:
$$
x_1=x_0-\frac{f(x_0)}{f’(x_0)}
$$
这样就得到根的新的近似值。进一步迭代:
$$
x_{k+1}=x_k-\frac{f(x_k)}{f’(x_k)}
$$
当前后两步解的差值达到某个精度时,就认为该解就是方程的根。

多变量求根

回到有限元法中,因为场函数由分布在节点上的有限元函数值所代替,所以此时变量就有多个。那么,对应于单变量情形,可得出多变量时的表达式:
$$
\mathbf{J}(\vec{u}_0) (\vec{u}_1-\vec{u}_0) = -\vec{R}(\vec{u}_0)
$$
其中,$\mathbf{J}(\vec{u}_0)$是当前迭代步的雅各比矩阵。上方带箭头的变量表示其是一个矢量,矢量大小就是节点的个数。
通用的迭代格式就是:
$$
\mathbf{J}(\vec{u}_k) (\vec{u}_{k+1}-\vec{u}_k) = -\vec{R}(\vec{u}_k)
$$
其中雅各比矩阵的具体形式为:
$$
J_{ij}(\vec{u}_k) = \frac{\partial R_i(\vec{u}_k)}{\partial u_j}
$$
在求雅各比矩阵时,有两个基本公式是非常重要的:
$$
\frac{\partial u_h}{\partial u_j} =
\sum_k\frac{\partial }{\partial u_j}(u_k \phi_k) = \phi_j
\qquad
\frac{\partial (\nabla u_h)}{\partial u_j} =
\sum_k \frac{\partial }{\partial u_j}(u_k \nabla \phi_k)=\nabla \phi_j
$$

牛刀小试

对于一个经典的对流-扩散方程:
$$
-\nabla\cdot k\nabla u + \vec{\beta} \cdot \nabla u = f
$$
其残差向量的第i个分量为:
$$
R_i(u_h) = (\nabla\psi_i, k\nabla u_h)-\langle\psi_i, k\nabla u_h\cdot \hat{n} \rangle + (\psi_i, \vec{\beta} \cdot \nabla u_h)-(\psi_i, f)
$$
雅各比矩阵的某个元素为:
$$
J_{ij}(u_h)= (\nabla\psi_i,\frac{\partial k}{\partial u_j}\nabla u_h) +(\nabla\psi_i, k \nabla \phi_j ) - \langle\psi_i, \frac{\partial k}{\partial u_j}\nabla u_h\cdot \hat{n} \rangle - \langle\psi_i, k\nabla \phi_j\cdot \hat{n} \rangle + (\psi_i, \frac{\partial \vec{\beta}}{\partial u_j} \cdot\nabla u_h) + (\psi_i, \vec{\beta} \cdot \nabla \phi_j) - (\psi_i, \frac{\partial f}{\partial u_j})
$$
注意,这里假定扩散系数、对流速度、源项等都是变量,所以形式会比较复杂。
尤其是对于多个方程相互耦合,或材料属性比较复杂的情形,雅各比矩阵的计算会很困难。

链式法则

有时候,可以应用链式法则来简化一下,比如:
$$
\frac{\partial f}{\partial u_j} = \frac{\partial f}{\partial u_h} \frac{\partial u_h}{\partial u_j}=\frac{\partial f}{\partial u_h} \phi_j
$$
如果f的表达式已知,比如$f(u) = \sin(u)$,那么:
$$
\frac{\partial f}{\partial u_j} = \cos(u_h) \phi_j
$$

JFNK算法

JFNK算法,全称是Jacobian Free Newton Krylov methods,是数值求解非线性问题的一种先进方法,其核心思想是将Newton非线性迭代法嵌入到Krylov空间法求解线性代数方程组的过程中,其显著优点是避免了传
统Newton迭代法中的Jacobian矩阵生成环节,有利于降低内存占用率,缩短计算时长。
JFNK算法将Newton迭代法与Krylov空间法结合的方式是将Newton迭代法中的Jacobian矩阵$\mathbf{J}$与Krylov空间法的解向量$\vec{v}$之间进行向量积操作,其近似为:
$$
\mathbf{J}\vec{v} \approx \frac{\vec{R}(\vec{u} + \epsilon\vec{v}) - \vec{R}(\vec{u})}{\epsilon}
$$
这个算法的优点有:

  • 无需计算偏导数来得到雅各比矩阵
  • 无需直接计算雅各比矩阵
  • 无需空间来存储雅各比矩阵

调幅分解所用的Cahn-Hilliard方程

这里将含有四阶导数的原CH方程,拆分成两个,这样每个都只含有二阶导数,易于求解,两个方程的变量分别是化学势和浓度。

化学势的残差

其表达式及其求解内核分为两个:

第一项

$$
(\frac{\partial c}{\partial t},\psi)
$$
该项中变量是化学势,耦合的变量是浓度,使用的内核是CoupledTimeDerivative。

1
2
3
4
5
6
7
8
9
10
11
Real
CoupledTimeDerivative::computeQpResidual()
{
return _test[_i][_qp] * _v_dot[_qp];
}

Real
CoupledTimeDerivative::computeQpJacobian()
{
return 0.0;
}

第二项

$$
(M\nabla\mu,\nabla \psi)
$$
该项中变量是化学势,使用的内核是SplitCHWRes,实际使用的是SplitCHWResBase:

1
2
3
4
5
6
7
8
9
10
11
12
13
template<typename T>
Real
SplitCHWResBase<T>::computeQpResidual()
{
return _mob[_qp] * _grad_u[_qp] * _grad_test[_i][_qp];
}

template<typename T>
Real
SplitCHWResBase<T>::computeQpJacobian()
{
return _mob[_qp] * _grad_phi[_j][_qp] * _grad_test[_i][_qp];
}

浓度的残差

残差表达式为:
$$
(\nabla c, \nabla(\kappa \psi))+((\frac{\partial f_{loc}}{\partial c}+\frac{\partial E_d}{\partial c} - \mu ), \psi)
$$
变量是浓度,还需耦合化学势$\mu$,使用的内核是SplitCHParsed,实际使用的内核是SplitCHCRes:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
 Real
SplitCHCRes::computeQpResidual()
{
Real residual = SplitCHBase::computeQpResidual(); //(f_prime_zero+e_prime)*_test[_i][_qp] from SplitCHBase

residual += -_w[_qp] * _test[_i][_qp];
residual += _kappa[_qp] * _grad_u[_qp] * _grad_test[_i][_qp];

return residual;
}

Real
SplitCHCRes::computeQpJacobian()
{
Real jacobian = SplitCHBase::computeQpJacobian(); //(df_prime_zero_dc+de_prime_dc)*_test[_i][_qp]; from SplitCHBase

jacobian += _kappa[_qp] * _grad_phi[_j][_qp] * _grad_test[_i][_qp];

return jacobian;
}
```

注意,在computeQpResidual函数中,可以很容易地找出第一项和第三项的计算过程,但对于第二项,可能不容易发现,其实第二项的计算放在了最前面:

```cpp
Real residual = SplitCHBase::computeQpResidual(); //(f_prime_zero+e_prime)*_test[_i][_qp] from SplitCHBase

然后:

1
2
3
4
5
6
7
8
9
10
Real
SplitCHBase::computeQpResidual()
{
Real f_prime_zero = computeDFDC(Residual);
Real e_prime = computeDEDC(Residual);

Real residual = (f_prime_zero + e_prime) *_test[_i][_qp];

return residual;
}

然后computeDFDC和computeDEDC就是计算自由能密度和其他能量对浓度的一阶导数,即:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
Real
SplitCHParsed::computeDFDC(PFFunctionType type)
{
switch (type)
{
case Residual:
return _dFdc[_qp];

case Jacobian:
return _d2Fdc2[_qp] * _phi[_j][_qp];
}

mooseError("Internal error");
}

实际上,

1
2
_dFdc(getMaterialPropertyDerivative<Real>("f_name", _var.name())),
_d2Fdc2(getMaterialPropertyDerivative<Real>("f_name", _var.name(), _var.name()))

注意取一阶导数和二阶导数时,是由后面的参数个数来控制。

具体的自由能形式则是在输入文件中输入:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
[Materials]
[./local_energy]
# Defines the function for the local free energy density as given in the
# problem, then converts units and adds scaling factor.
type = DerivativeParsedMaterial
block = 0
f_name = f_loc
args = c
constant_names = 'A B C D E F G eV_J d'
constant_expressions = '-2.446831e+04 -2.827533e+04 4.167994e+03 7.052907e+03
1.208993e+04 2.568625e+03 -2.354293e+03
6.24150934e+18 1e-27'
function = 'eV_J*d*(A*c+B*(1-c)+C*c*log(c)+D*(1-c)*log(1-c)+
E*c*(1-c)+F*c*(1-c)*(2*c-1)+G*c*(1-c)*(2*c-1)^2)'
[../]
[]

枝晶生长所用的Kobayashi模型

相场的残差

包括以下几项:

第一项

$$
(\frac{\partial w}{\partial t},\psi)
$$
内核使用TimeDerivative:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
#include "TimeDerivative.h"
#include "Assembly.h"

// libmesh includes
#include "libmesh/quadrature.h"

template<>
InputParameters validParams<TimeDerivative>()
{
InputParameters params = validParams<TimeKernel>();
params.addParam<bool>("lumping", false, "True for mass matrix lumping, false otherwise");
return params;
}

TimeDerivative::TimeDerivative(const InputParameters & parameters) :
TimeKernel(parameters),
_lumping(getParam<bool>("lumping"))
{
}

Real
TimeDerivative::computeQpResidual()
{
return _test[_i][_qp]*_u_dot[_qp];
}

Real
TimeDerivative::computeQpJacobian()
{
return _test[_i][_qp]*_phi[_j][_qp]*_du_dot_du[_qp];
}

void
TimeDerivative::computeJacobian()
{
if (_lumping)
{
DenseMatrix<Number> & ke = _assembly.jacobianBlock(_var.number(), _var.number());

for (_i = 0; _i < _test.size(); _i++)
for (_j = 0; _j < _phi.size(); _j++)
for (_qp = 0; _qp < _qrule->n_points(); _qp++)
ke(_i, _i) += _JxW[_qp] * _coord[_qp] * computeQpJacobian();
}
else
TimeKernel::computeJacobian();
}

第二项

注意将残差都移动到方程左端,注意符号变化。
$$
(-L\epsilon\epsilon\prime\frac{\partial w}{\partial y},\frac{\partial \psi}{\partial x})+(L\epsilon\epsilon\prime\frac{\partial w}{\partial x},\frac{\partial \psi}{\partial y})
$$
这两项使用内核ACInterfaceKobayashi1.h,其实际使用KernelGrad:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
RealGradient
ACInterfaceKobayashi1::precomputeQpResidual()
{
// Set modified gradient vector
const RealGradient v(- _grad_u[_qp](1), _grad_u[_qp](0), 0);

// Define anisotropic interface residual
return _eps[_qp] * _deps[_qp] * _L[_qp] * v;
}
...
void
KernelGrad::computeResidual()
{
DenseVector<Number> & re = _assembly.residualBlock(_var.number());
_local_re.resize(re.size());
_local_re.zero();

const unsigned int n_test = _test.size();
for (_qp = 0; _qp < _qrule->n_points(); _qp++)
{
RealGradient value = precomputeQpResidual() * _JxW[_qp] * _coord[_qp];
for (_i = 0; _i < n_test; _i++) // target for auto vectorization
_local_re(_i) += value * _grad_test[_i][_qp];
}

re += _local_re;

if (_has_save_in)
{
Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
for (const auto & var : _save_in)
var->sys().solution().add_vector(_local_re, var->dofIndices());
}
}

第三项

$$
(L\epsilon^2\nabla w,\nabla \psi)
$$
注意,内核中不是直接使用该形式,该形式是自由能对相场参量的一阶导数,实际程序中使用的也是一阶导数,所以需要输入的是自由能,即该形式的积分。
其使用的内核是ACInterfaceKobayashi2,实际也是KernelGrad:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
RealGradient
ACInterfaceKobayashi2::precomputeQpResidual()
{
// Set interfacial part of residual
return _eps[_qp] * _eps[_qp] * _L[_qp] * _grad_u[_qp];
}
...
KernelGrad::computeResidual()
{
DenseVector<Number> & re = _assembly.residualBlock(_var.number());
_local_re.resize(re.size());
_local_re.zero();

const unsigned int n_test = _test.size();
for (_qp = 0; _qp < _qrule->n_points(); _qp++)
{
RealGradient value = precomputeQpResidual() * _JxW[_qp] * _coord[_qp];
for (_i = 0; _i < n_test; _i++) // target for auto vectorization
_local_re(_i) += value * _grad_test[_i][_qp];
}

re += _local_re;

if (_has_save_in)
{
Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
for (const auto & var : _save_in)
var->sys().solution().add_vector(_local_re, var->dofIndices());
}
}

(4)
$$
(-Lw(1-w)(w-0.5+m),\psi)
$$
使用的内核是AllenCahn,实际间接使用了ACBulk,最终使用KernelValue:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
...
_dFdEta(getMaterialPropertyDerivative<Real>("f_name", _var.name())),
_d2FdEta2(getMaterialPropertyDerivative<Real>("f_name", _var.name(), _var.name())),
...
Real
AllenCahn::computeDFDOP(PFFunctionType type)
{
switch (type)
{
case Residual:
return _dFdEta[_qp];

case Jacobian:
return _d2FdEta2[_qp] * _phi[_j][_qp];
}

mooseError("Internal error");
}
...
template<typename T>
Real
ACBulk<T>::precomputeQpResidual()
{
// Get free energy derivative from function
Real dFdop = computeDFDOP(Residual);

// Set residual
return _L[_qp] * dFdop;
}
...
void
KernelValue::computeResidual()
{
DenseVector<Number> & re = _assembly.residualBlock(_var.number());
_local_re.resize(re.size());
_local_re.zero();

const unsigned int n_test = _test.size();
for (_qp = 0; _qp < _qrule->n_points(); _qp++)
{
Real value = precomputeQpResidual() * _JxW[_qp] * _coord[_qp];
for (_i = 0; _i < n_test; _i++) // target for auto vectorization
_local_re(_i) += value * _test[_i][_qp];
}

re += _local_re;

if (_has_save_in)
{
Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
for (const auto & var : _save_in)
var->sys().solution().add_vector(_local_re, var->dofIndices());
}
}

温度场的残差

由以下几项构成:

第一项

$$
(\frac{\partial T}{\partial t},\psi)
$$
所以使用的内核是TimeDerivative。

第二项

$$
(\nabla T,\nabla \psi)
$$
所以使用的内核是Diffusion。

第三项

$$
(-K\frac{\partial w}{\partial t},\psi)
$$
使用的内核是CoefCoupledTimeDerivative,实际就是CoupledTimeDerivative再乘以一个系数:

1
2
3
4
5
6
7
8
9
10
Real
CoefCoupledTimeDerivative::computeQpResidual()
{
return CoupledTimeDerivative::computeQpResidual() * _coef;
}
Real
CoupledTimeDerivative::computeQpResidual()
{
return _test[_i][_qp] * _v_dot[_qp];
}