Palabos源码：momentTemplates3D

最新推荐文章于 2021-04-10 19:15:13 发布

方雨岚(Yulan Fang)

最新推荐文章于 2021-04-10 19:15:13 发布

阅读量233

点赞数

分类专栏： Palabos进阶学习

本文链接：https://blog.csdn.net/qq_40259141/article/details/111568421

版权

Palabos进阶学习专栏收录该内容

37 篇文章 47 订阅

订阅专栏

这是比较底层的代码，基本上我们在palabos的各种操作都是基于数据处理器上的，而数据处理器则是通过各种方式调用这些底层的代码进行基于分布函数f的运算。
以D3Q19为例，介绍一些常见的功能。

typedef descriptors::D3Q19DescriptorBase Descriptor;

get_rhoBar

T rhoBar = f[0] + f[1] + f[2] + f[3] + f[4]
                   + f[5] + f[6] + f[7] + f[8]
                   + f[9] + f[10] + f[11] + f[12]
                   + f[13] + f[14] + f[15] + f[16]
                   + f[17] + f[18];
return rhoBar;

get_j

T surfX_M1, surfX_P1, surfY_M1, surfY_P1, surfZ_M1, surfZ_P1;

surfX_M1 = f[1] + f[4] + f[5] + f[6] + f[7];
surfX_P1 = f[10] + f[13] + f[14] + f[15] + f[16];

surfY_M1 = f[2] + f[4] + f[8] + f[9] + f[14];
surfY_P1 = f[5] + f[11] + f[13] + f[17] + f[18];

surfZ_M1 = f[3] + f[6] + f[8] + f[16] + f[18];
surfZ_P1 = f[7] + f[9] + f[12] + f[15] + f[17];

j[0]  = ( surfX_P1 - surfX_M1 );
j[1]  = ( surfY_P1 - surfY_M1 );
j[2]  = ( surfZ_P1 - surfZ_M1 );

get_rhoBar_j

T surfX_M1, surfX_0, surfX_P1,
  surfY_M1, surfY_P1, surfZ_M1, surfZ_P1;

partial_rho(f, surfX_M1, surfX_0, surfX_P1,
                  surfY_M1, surfY_P1, surfZ_M1, surfZ_P1);

rhoBar = surfX_M1 + surfX_0 + surfX_P1;

j[0]  = ( surfX_P1 - surfX_M1 );
j[1]  = ( surfY_P1 - surfY_M1 );
j[2]  = ( surfZ_P1 - surfZ_M1 );

partial_rho

static void partial_rho ( Array<T,Descriptor::q> const& f,
                          T& surfX_M1, T& surfX_P1,
                          T& surfY_M1, T& surfY_P1, T& surfZ_M1, T& surfZ_P1 )
{
    surfX_M1 = f[1] + f[4] + f[5] + f[6] + f[7];
    surfX_P1 = f[10] + f[13] + f[14] + f[15] + f[16];

    surfY_M1 = f[2] + f[4] + f[8] + f[9] + f[14];
    surfY_P1 = f[5] + f[11] + f[13] + f[17] + f[18];

    surfZ_M1 = f[3] + f[6] + f[8] + f[16] + f[18];
    surfZ_P1 = f[7] + f[9] + f[12] + f[15] + f[17];

}

static void partial_rho ( Array<T,Descriptor::q> const& f,
                          T& surfX_M1, T& surfX_0, T& surfX_P1,
                          T& surfY_M1, T& surfY_P1, T& surfZ_M1, T& surfZ_P1 )
{
    partial_rho(f, surfX_M1, surfX_P1, surfY_M1, surfY_P1, surfZ_M1, surfZ_P1);
    surfX_0  = f[0] + f[2] + f[3] + f[8] +
               f[9] + f[11] + f[12] + f[17] + f[18];
}

compute_uLb

get_j(f, uLb);
T invRho = Descriptor::invRho(get_rhoBar(f));
uLb[0] *= invRho;
uLb[1] *= invRho;
uLb[2] *= invRho;

compute_rho_uLb

T rhoBar;
get_rhoBar_j(f, rhoBar, uLb);
T invRho = Descriptor::invRho(rhoBar);
rho = Descriptor::fullRho(rhoBar);
uLb[0] *= invRho;
uLb[1] *= invRho;
uLb[2] *= invRho;

compute_rhoBar_j_PiNeq

static void compute_rhoBar_j_PiNeq(Array<T,Descriptor::q> const& f, T& rhoBar, Array<T,3>& j, Array<T,6>& PiNeq)
{
    typedef SymmetricTensorImpl<T,3> S;

    T surfX_M1, surfX_0, surfX_P1,
      surfY_M1, surfY_P1, surfZ_M1, surfZ_P1;

    partial_rho(f, surfX_M1, surfX_0, surfX_P1,
                   surfY_M1, surfY_P1, surfZ_M1, surfZ_P1);
    rhoBar = surfX_M1 + surfX_0 + surfX_P1;
    j[0]  = ( surfX_P1 - surfX_M1 );
    j[1]  = ( surfY_P1 - surfY_M1 );
    j[2]  = ( surfZ_P1 - surfZ_M1 );
    T invRho = Descriptor::invRho(rhoBar);

    PiNeq[S::xx] = surfX_P1+surfX_M1 - Descriptor::cs2*rhoBar - invRho*j[0]*j[0];
    PiNeq[S::yy] = surfY_P1+surfY_M1 - Descriptor::cs2*rhoBar - invRho*j[1]*j[1];
    PiNeq[S::zz] = surfZ_P1+surfZ_M1 - Descriptor::cs2*rhoBar - invRho*j[2]*j[2];

    PiNeq[S::xy] = f[4] - f[5] + f[13] - f[14] - invRho*j[0]*j[1];
    PiNeq[S::xz] = f[6] - f[7] + f[15] - f[16] - invRho*j[0]*j[2];
    PiNeq[S::yz] = f[8] - f[9] + f[17] - f[18] - invRho*j[1]*j[2];
}

static void compute_rhoBar_j_PiNeq(Array<T,Descriptor::q> const& f, T& rhoBar, Array<T,3>& j, Array<T,6>& PiNeq, T invRho)
{
    typedef SymmetricTensorImpl<T,3> S;

    T surfX_M1, surfX_0, surfX_P1,
      surfY_M1, surfY_P1, surfZ_M1, surfZ_P1;

    partial_rho(f, surfX_M1, surfX_0, surfX_P1,
                   surfY_M1, surfY_P1, surfZ_M1, surfZ_P1);
    rhoBar = surfX_M1 + surfX_0 + surfX_P1;
    j[0]  = ( surfX_P1 - surfX_M1 );
    j[1]  = ( surfY_P1 - surfY_M1 );
    j[2]  = ( surfZ_P1 - surfZ_M1 );

    PiNeq[S::xx] = surfX_P1+surfX_M1 - Descriptor::cs2*rhoBar - invRho*j[0]*j[0];
    PiNeq[S::yy] = surfY_P1+surfY_M1 - Descriptor::cs2*rhoBar - invRho*j[1]*j[1];
    PiNeq[S::zz] = surfZ_P1+surfZ_M1 - Descriptor::cs2*rhoBar - invRho*j[2]*j[2];

    PiNeq[S::xy] = f[4] - f[5] + f[13] - f[14] - invRho*j[0]*j[1];
    PiNeq[S::xz] = f[6] - f[7] + f[15] - f[16] - invRho*j[0]*j[2];
    PiNeq[S::yz] = f[8] - f[9] + f[17] - f[18] - invRho*j[1]*j[2];
}

compute_P

typedef SymmetricTensorImpl<T,3> S;

T invRho = Descriptor::invRho(rhoBar);
T surfX_M1, surfX_P1, surfY_M1, surfY_P1, surfZ_M1, surfZ_P1;
partial_rho(f, surfX_M1, surfX_P1, surfY_M1, surfY_P1, surfZ_M1, surfZ_P1);

P[S::xx] = surfX_P1+surfX_M1 - invRho*j[0]*j[0];
P[S::yy] = surfY_P1+surfY_M1 - invRho*j[1]*j[1];
P[S::zz] = surfZ_P1+surfZ_M1 - invRho*j[2]*j[2];

P[S::xy] = f[4] - f[5] + f[13] - f[14] - invRho*j[0]*j[1];
P[S::xz] = f[6] - f[7] + f[15] - f[16] - invRho*j[0]*j[2];
P[S::yz] = f[8] - f[9] + f[17] - f[18] - invRho*j[1]*j[2];

modifyJ

T rhoBar, oldJ[3];
get_rhoBar_j(f, rhoBar, oldJ);
const T oldJSqr = VectorTemplateImpl<T,3>::normSqr(oldJ);
const T newJSqr = VectorTemplateImpl<T,3>::normSqr(newJ);
for (plint iPop=0; iPop<Descriptor::q; ++iPop) {
    f[iPop] = f[iPop]
                     - equilibrium(iPop, rhoBar, oldJ, oldJSqr)
                     + equilibrium(iPop, rhoBar, newJ, newJSqr);