openMP交替方向隐式方法

Question

嗨，我正在尝试用openMP来实现这段代码的计算。采用有限差分隐式方法计算流体动力学涡量。我使用交替方向隐式方法来做到这一点。

我想加快它的执行速度。(此处为Nx=Ny=100)

问题是，以这种方式使用openMp会降低代码的速度，而不是加快代码的速度。我已经尝试指定共享变量，但这没有多大帮助。有什么想法吗？

万事如意

void ADI(double vort[][Ny], double psi[][Ny], double n[][Ny], 
double cls[][Ny],double AAx[], double BBx[], double CCx[], double DDx[],
double AAy[], double BBy[], double CCy[], double DDy[],
double cx[][Ny], double cy[][Ny], double epsx[][Ny], double epsy[][Ny],
double vortx[], double vorty[Ny-2], double dx, double Dxs, double coefMass, 
double coefMasCls)
{
    ////////////calcul sur y////////////
    //calcul coef ADI
    int i=0, j=0;

    #pragma omp parallel for private(Dxs,i) shared(psi,vort)
    for (i=0; i<Nx; i++) //Boundary condition sur x
    {
        vort[i][0]=(psi[i][0]-psi[i][1])*2/Dxs;
        vort[i][Ny-1] = (psi[i][Ny-1]-psi[i][Ny-2])*2/Dxs; 
    }

    #pragma omp parallel for private(Dxs,j) shared(psi,vort)
    for (j=0; j<Ny; j++) //Boundary condition
    {
        vort[0][j] = (psi[0][j]-psi[1][j])*2/Dxs;
        vort[Nx-1][j] = (psi[Nx-1][j]-psi[Nx-2][j])*2/Dxs;

    }

    for (j=1; j<Ny-1; j++) //interior points
    {
        #pragma omp parallel for private(coefMasCls,coefMasCls,i) shared(psi,vort,n,cls)
        for (i=1; i<Nx-1; i++) //interior points
        {

            vort[i][j] = vort[i][j] - coefMass * (n[i+1][j]-n[i-1][j])- coefMasCls * (cls[i+1][j]-cls[i-1][j]);;

        }
        //i=0;
        //vort[i][j] = vort[i][j] + coefMass*(n[1][j]-n[1][j]);
        //i=Nx-1;
        //vort[i][j] = vort[i][j] + coefMass*(n[Nx-2][j]-n[Nx-2][j]);

    }


    for (i=1; i<Nx-1; i++) //interior points
    {
        for (j=1; j<Ny-1; j++) //interior points
        {


            AAy[j] = -.5 * ( .5 * (1 + epsy[i][j]) * cy[i][j-1] + dx);

            BBy[j] = 1 + dx + .5 * epsy[i][j] * cy[i][j];

            CCy[j] = .5 * ( .5 * ( 1 - epsy[i][j] ) * cy[i][j+1] - dx);

            DDy[j] = .5 * (.5 * ( 1 + epsx[i][j] ) * cx[i-1][j] + dx ) * vort[i-1][j]
            + ( 1 - dx - .5 * epsx[i][j] * cx[i][j] ) * vort[i][j]
            + .5 * (- .5 * ( 1 - epsx[i][j] ) * cx[i+1][j] + dx ) * vort[i+1][j];


            vorty[j] = vort[i][j];
        }


        DDy[1]=DDy[1] - AAy[1] * vort[i][0];  //the AA[0] are not taken into account in the tridiag methode. Include it in the second hand
        DDy[Ny-2]=DDy[Ny-2] - CCy[Ny-2]* vort[i][Ny-1]; //moving boundary condition
        //DDy[Ny-3]= DDy[Ny-3]; //vorticity nul on the free slip boundary condition


        tridiag(AAy, BBy, CCy, DDy, vorty, Ny-1); //ne calcul pas le point en 0 et en Ny-1

        for (j=1; j<Ny-1; j++)
        {
            vort[i][j]=vorty[j];
        }

    }

    ////////////calcul sur x //////////
    //calcul coef ADI

    for (j=1; j<Ny-1; j++)
    {

        for (i=1; i<Nx-1; i++)

        {

            AAx[i] = -.5* ( .5 * ( 1 + epsx[i][j] ) * cx[i-1][j] + dx );
            BBx[i] = 1 + dx + .5 * epsx[i][j] * cx[i][j];
            CCx[i] = .5 * ( .5 * ( 1 - epsx[i][j] ) * cx[i+1][j] - dx) ;

            DDx[i]= .5 * ( .5 * ( 1 + epsy[i][j] ) * cy[i][j-1] + dx ) * vort[i][j-1]
            + ( 1 - dx - .5 * epsy[i][j] * cy[i][j] ) * vort[i][j]
            + .5 * (-.5 * ( 1 - epsy[i][j] ) * cy[i][j+1] + dx ) * vort[i][j+1];

            vortx[i]=vort[i][j];

        }


        DDx[1] = DDx[1] - AAx[1]* vort[0][j];
        DDx[Nx-2] = DDx[Nx-2] - CCx[Nx-2] * vort[Nx-1][j];

        tridiag(AAx, BBx, CCx, DDx, vortx, Nx-1); //ne calcul pas le point en 0 et en Nx-1

        for (i=1; i<Nx-1; i++)
        {
            vort[i][j]=vortx[i];

        }


    }

}

Answer 1

首先要做的确实是隔离哪些循环并行化有最坏的影响，但最后一个循环看起来很像是在经历缓存颠簸。稍微简化一下结构：

double vort[Nx][Ny];
// ...
for (int j=1; j<Ny-1; ++j) {
    #pragma omp parallel for
    for (int i=1; i<Nx-1; ++i) {
        vort[i][j] -= f(i, j);
    }
}

任何给定的线程都将在偏移量j+k*Ny，j+(k+1)*Ny，j+(k+2)*Ny等处依次读取和更新vort中的值，这取决于for循环是如何跨线程分块的。这些访问中的每一次都会引入一个缓存线的数据来更新8个字节。当外部循环再次启动时，您刚刚访问的数据很可能都不会仍在缓存中。

在所有条件相同的情况下，如果你能安排你的数组访问，使你在最小步长的方向上移动(对于C数组，这是最后一个索引)，你的缓存行为将会更好。对于维度大小为100的数组，数组可能不会太大，这会产生巨大的差异。例如，对于Nx，Ny = 1000，以“错误的方式”访问数组可能是灾难性的。

这会使串行代码的性能变差，但我认为添加线程只会使情况变得更糟。

尽管如此，在每个内部循环中完成的计算量都非常小；无论如何，您都很有可能会受到内存带宽的限制。

附录

明确地说，“正确”的循环访问将如下所示：

for (int i=1; i<Nx-1; ++i) {
    for (int j=1; j<Ny-1; ++j) {
        vort[i][j] -= f(i, j);
    }
}

要使其并行化，可以通过使用collapse指令让编译器更好地跨线程分块数据：

#pragma omp parallel for collapse(2)
for (int i=1; i<Nx-1; ++i) {
    for (int j=1; j<Ny-1; ++j) {
        vort[i][j] -= f(i, j);
    }
}

最后，为了避免假共享(线程相互踩在对方的缓存线上)，最好确保数组的两个相邻行不共享同一缓存线中的数据。可以确保每一行都与内存中缓存行大小的倍数对齐，或者更简单地将填充添加到每一行的末尾。

double vort[Nx][Ny+8]; // 8 doubles ~ 64 bytes

(假设缓存线为64字节，这应该足够了。)