交错网格上的几何多重网格方法

发表于 编辑

对于规则网格而言，未知量有多种不同的布置方式，以此为依据又可分为顶点中心网格（vertex-centered grid）、单元中心网格（cell-centered grid）和交错网格（staggered grid）。使用多重网格方法求解问题时，不同类型的网格需要采用不同的限制算子(restriction operator)和延拓算子(prolongation(interpolation) operator)。限制算子和延拓算子统称为传递算子(transfer operator)。

对于顶点中心网格，限制算子通常使用完全加权算子(full weighting operator)，此外还有半加权算子(half weighting operator)，injection operator，half injection operator等；对于单元中心网格，限制算子通常使用平均法，二维使用四点平均法，三维使用八点平均法。

由于非线性项的存在，在求解稳态不可压Navier-Stokes方程时采用了FAS（Full Approximation Scheme）。因此，需要将未知量从细网格限制到粗网格上。对于限制算子，使用邻近网格点的未知数的平均值来定义粗网格上的近似值。

uˉH=12[1.1]uˉh,vˉH=12[1⋅]vˉh,pˉH=14[11⋅11]pˉh,\begin{aligned} & \bar{u}_H=\frac{1}{2}\left[\begin{array}{l} 1 \\ . \\ 1 \end{array}\right] \bar{u}_h, \\ & \bar{v}_H=\frac{1}{2}\left[\begin{array}{ll} 1 & \cdot \end{array}\right] \bar{v}_h, \\ & \bar{p}_H=\frac{1}{4}\left[\begin{array}{lll} 1 & & 1 \\ & \cdot & \\ 1 & & 1 \end{array}\right] \bar{p}_h, \\ & \end{aligned}

从粗网格到细网格上残差的限制算子可以表示成

dHu=18[121121]dhu,dHv=18[112⋅211]dhv,dHp=14[1111]dhp.\begin{aligned} & d_H^u=\frac{1}{8}\left[\begin{array}{lll} 1 & 2 & 1 \\ 1 & 2 & 1 \end{array}\right] d_h^u, \\ & d_H^v=\frac{1}{8}\left[\begin{array}{lll} 1 & & 1 \\ 2 & \cdot & 2 \\ 1 & & 1 \end{array}\right] d_h^v, \\ & d_H^p=\frac{1}{4}\left[\begin{array}{lll} 1 & 1 \\ 1 & & 1 \end{array}\right] d_h^p . \end{aligned}

箱形松弛方法(box relaxation method): 对于任意网格单元(i,j)(i,j)，存在五个未知量：ui−1,j−1/2u_{i-1,j-1/2}、vi−1/2,j−1v_{i-1/2,j-1}、pi−1/2,j−1/2p_{i-1/2,j-1/2}、ui,j−1/2u_{i,j-1/2}和vi−1/2,jv_{i-1/2,j}，箱形松弛方法将这五个变量视为一个集体，利用高斯-塞德尔迭代思想，使用其他未知量来更新这五个未知量。一次完整的迭代需要遍历所有网格单元，更新每个单元上的未知量。

xn+1−xnΔt=un\frac{\mathbf{x}^{n+1}-\mathbf{x}^n}{\Delta t}=\mathbf{u}^{n}

udn=u(xn,t)\mathbf{u}^n_d=\mathbf{u}(\mathbf{x}^n,t)

线性插值算子

延拓算子可使用线性插值算子。线性插值推广到二维、三维，可以得到双线性插值和三线性插值，它们的公式分别为：

f(x,y,z)=∑i=01∑j=01∑k=01fi,j,kLi(x)Lj(y)Lk(z),f(x,y,z)=∑i=01∑j=01fi,j,kLi(x)Lj(y),f(x,y,z)=∑i=01fi,j,kLi(x).\begin{aligned} & f(x, y, z)=\sum_{i=0}^1 \sum_{j=0}^1 \sum_{k=0}^1 f_{i, j, k} L_i(x) L_j(y) L_k(z), \\ & f(x, y, z)=\sum_{i=0}^1 \sum_{j=0}^1 f_{i, j, k} L_i(x) L_j(y),\\ & f(x, y, z)=\sum_{i=0}^1 f_{i, j, k} L_i(x). \end{aligned}

其中，fi,j,kf_{i,j,k}为函数f(x,y,z)f(x,y,z)在线段、矩形或长方体顶点处的值，x,y,zx,y,z为区域内的坐标，L0L_0和L1L_1为

L0(x)=1.0−x,L1(x)=x.L_0(x)=1.0-x,\\ L_1(x)=x.

因此只要知道密网格和细网格的位置相对关系，就能确定延拓算子，三个分量的延拓算子可以定义成
速度分量u的插值
速度分量v的插值
压强p的插值

可通过如下C++代码来检验。

#include <array>
#include <functional>
#include <iostream>

// 考虑GPU设备运行，使用函数指针而不是lamgbda函数
double L0(double x) { return 1.0 - x; }
double L1(double x) { return x; }
double (*L[2])(double) = {L0, L1};

int test_linears()
{
    double f[2];
    // f(x) = x+1
    f[0] = 0; // (0)
    f[1] = 1; // (1)

    // std::array<std::function<double(double)>, 2> L = {{
    //     [](double x) { return (1.0-x); },
    //     [](double x) { return (x); }
    // }};

    double sum = 0.0;
    double x = 0.1;
    double y = 0.3;
    for (int i = 0; i < 2; i++)
    {
        sum += f[i] * L[i](x);
        std::cout << sum << std::endl;
    }
    return 0;
}

int test_bilinears()
{
    double f[2][2];
    // f(x,y) = x+y+2xy
    f[0][0] = 0; // (0,0)
    f[1][0] = 1; // (1,0)
    f[0][1] = 1; // (0,1)
    f[1][1] = 4; // (1,1)

    // std::array<std::function<double(double)>, 2> L = {{
    //     [](double x) { return (1.0-x); },
    //     [](double x) { return (x); }
    // }};

    double sum = 0.0;
    double x = 0.1;
    double y = 0.3;
    for (int i = 0; i < 2; i++)
        for (int j = 0; j < 2; j++)
        {
            sum += f[i][j] * L[i](x) * L[j](y);
            std::cout << sum << std::endl;
        }
    return 0;
}

// 三线性插值
double f_(double x, double y, double z)
{
    return x + y + z + x * y + x * z + y * z + x * y * z + 1;
}
int main()
{
    double f[2][2][2];
    // f(x,y) = x+y+z+xy+xz+yz+xyz+1
    f[0][0][0] = f_(0, 0, 0); //
    f[1][0][0] = f_(1, 0, 0); //
    f[0][1][0] = f_(0, 1, 0); //
    f[1][1][0] = f_(1, 1, 0); //
    f[0][0][1] = f_(0, 0, 1); //
    f[1][0][1] = f_(1, 0, 1); //
    f[0][1][1] = f_(0, 1, 1); //
    f[1][1][1] = f_(1, 1, 1); //

    // std::array<std::function<double(double)>, 2> L = {{
    //     [](double x) { return (1.0-x); },
    //     [](double x) { return (x); }
    // }};

    double sum = 0.0;
    double x = 0.2;
    double y = 0.3;
    double z = 0.3;
    for (int i = 0; i < 2; i++)
        for (int j = 0; j < 2; j++)
            for (int k = 0; k < 2; k++)
            {
                sum += f[i][j][k] * L[i](x) * L[j](y) * L[k](z);
            }

    std::cout << sum << std::endl;
    std::cout << f_(x, y, z) << std::endl;

    return 0;
}