用 OpenFOAM 进行张量计算

我们常说 OpenFOAM 是一个 C++ 类库，初学 OpenFOAM 时，我常对 C++ 类库这个概念感到困惑，OpenFOAM 不是一个用 C++ 写的 CFD 软件吗，怎么又成了 C++ 类库？其实 ，CFD 求解器只是 OpenFOAM 的一些高级应用。由于 OpenFOAM 包含了大量基础数学运算库，我们完全可以以一种更底层的方式使用 OpenFOAM，比如今天我们要谈到的张量计算。

在大学的第一堂 C 语言课上，我们便被告知，写程序之前，要先 include 头文件，比如这段 Hello World 代码。

#include <stdio.h>
int main()
{
   printf("Hello, World!");
   return 0;
}Code language: PHP (php)

这里的 stdio 就是标准输入输出库。而在 C++ 中，常用标准库 iostream 实现类似功能。

#include <iostream>
int main()
{
  std::cout << "Hello World!" << std::endl;
  return 0;
}
Code language: PHP (php)

再比如，如果要在 C++ 中使用数学函数，则要包含 cmath 库。OpenFOAM 提供了大量这样的数学运算库，比如用于张量计算的 tensor 和 symmTensor 等。

OpenFOAM 有一个现成的张量计算的例子，位置是 $FOAM_APP/test/tensor。我们拆开来说明下。

1. 相加

#include "tensor.H"
#include "IOstreams.H"
using namespace Foam;
int main()
{
    tensor t1(1, 2, 3, 4, 5, 6, 7, 8, 9);
    tensor t2(1, 2, 3, 1, 2, 3, 1, 2, 3);
    tensor t3 = t1 + t2;
    Info<< t3 << endl;
    return 0;
}
Code language: PHP (php)

t1、t2 和 t3 都是二阶张量，表示为三阶方阵。

2. 求逆

#include "tensor.H"
#include "IOstreams.H"
using namespace Foam;
int main()
{
    tensor t4(3,-2,1,-2,2,0,1, 0, 4);
    Info<< t4 << endl;
    Info<< inv(t4) << endl;
    return 0;
}
Code language: PHP (php)

inv() 是求逆（inverse）函数。

3. 特征值、特征向量与行列式

#include "tensor.H"
#include "IOstreams.H"
using namespace Foam;
int main()
{
    tensor t5(1,0,-4,0,5,4,-4,4,3);
    Info<< "tensor " << t5 << endl;
    vector e = eigenValues(t5);
    Info<< "eigenvalues " << e << endl;
    tensor ev = eigenVectors(t5);
    Info<< "eigenvectors " << ev << endl;
    Info<< "Check determinant " << e.x()*e.y()*e.z() << " " << det(t5) << endl;
    return 0;
}
Code language: PHP (php)

这里用到了一个数学性质——行列式等于特征值的积，详见 矩阵的特征：特征值，特征向量，行列式，trace。

此外，OpenFOAM 还提供了专门处理对称矩阵的库 symmTensor，详见 $FOAM_APP/test/tensor/Test_tensor.C。