c++ - Petsc 添加矩阵值

Question

我对 PETSC 有疑问。我在 matlab 中编写了一段代码，我正在尝试使用 PETSC 库将这段代码翻译成 C++。我正在为多相流编写流体动力学模拟，我正在尝试以一种简单的方式进行与此 matlab 操作等效的操作:

ut(i,j)=u(i,j)+(u(i+1,j)+u(i,j))^2-(u(i,j)+u(i-1,j))*(v(i, j) + v(i-1, j))

有没有一种方法可以做到这一点而不必调用 MATGETVALUES 函数 7 次？

Answer 1

如何将 u 和 v 存储为 vector ？您可能对 Petsc 的分布式阵列感兴趣，请参阅 DMDACreate2d()使用 2 个自由度 u 和 v。然后，函数 DMDAVecGetArrayDOF()将 DMDA 中的 vector 转换为数组。最后，你的行可以写成 ut[i][j][U]=ut[i][j][U]+.....+u[i-1][j][V ]).

下面这段代码基于this example显示如何完成:

static char help[] = "Just a piece of code. See ts/examples/tutorials/ex12.c\n";

#include <petscdm.h>
#include <petscdmda.h>

#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **argv)
{
    Vec                  x,r,xloc,rloc;                            
    PetscErrorCode       ierr;
    DM                   da;
    PetscScalar          ***xarray,***rarray;
    PetscInt             xm,ym,xs,ys,i,j,Mx,My;

    PetscInitialize(&argc,&argv,(char*)0,help);

    /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Create distributed array (DMDA) to manage parallel grid and vectors
  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    ierr = DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,-8,-8,PETSC_DECIDE,PETSC_DECIDE,
            2,1,NULL,NULL,&da);CHKERRQ(ierr);
    ierr = DMDASetFieldName(da,0,"u");CHKERRQ(ierr);
    ierr = DMDASetFieldName(da,1,"v");CHKERRQ(ierr);

    ierr = DMCreateGlobalVector(da,&x);CHKERRQ(ierr);
    ierr = VecDuplicate(x,&r);CHKERRQ(ierr);


    ierr = DMGetLocalVector(da,&xloc);CHKERRQ(ierr);
    ierr = DMGetLocalVector(da,&rloc);CHKERRQ(ierr);

    ierr = DMGlobalToLocalBegin(da,x,INSERT_VALUES,xloc);CHKERRQ(ierr);
    ierr = DMGlobalToLocalEnd(da,x,INSERT_VALUES,xloc);CHKERRQ(ierr);

    ierr = DMDAVecGetArrayDOF(da,xloc,&xarray);CHKERRQ(ierr);
    ierr = DMDAVecGetArrayDOF(da,rloc,&rarray);CHKERRQ(ierr);

    // Get global size
    ierr = DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);CHKERRQ(ierr);
    //Get local grid boundaries
    ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);
    for (j=ys; j<ys+ym; j++) {
        for (i=xs; i<xs+xm; i++) {

            if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
                //boundary conditions... Your job !
                continue;
            }
            // ut(i,j)=u(i,j)+(u(i+1,j)+u(i,j))^2-(u(i,j)+u(i-1,j))*(v(i, j) + v(i-1, j))
            rarray[j][i][0]=xarray[j][i][0]+pow(xarray[j][i+1][0]+xarray[j][i][0],2)-(xarray[j][i][0]+xarray[j][i-1][0])*(xarray[j][i][1]+xarray[j][i-1][1]);
            // v component
            rarray[j][i][1]=xarray[j][i][1]+pow(xarray[j+1][i][1]+xarray[j][i][1],2)-(xarray[j][i][1]+xarray[j-1][i][1])*(xarray[j][i][0]+xarray[j-1][i][0]);
        }
    }

    // the arrays are placed back in the vector.
    ierr = DMDAVecRestoreArrayDOF(da,xloc,&xarray);CHKERRQ(ierr);
    ierr = DMDAVecRestoreArrayDOF(da,rloc,&rarray);CHKERRQ(ierr);

    // local vector rloc is inserted in global vector.
    ierr = DMLocalToGlobalBegin(da,rloc,INSERT_VALUES,r);CHKERRQ(ierr);
    ierr = DMLocalToGlobalEnd(da,rloc,INSERT_VALUES,r);CHKERRQ(ierr);

    // local vector are destroyed.
    ierr = DMRestoreLocalVector(da,&xloc);CHKERRQ(ierr);
    ierr = DMRestoreLocalVector(da,&rloc);CHKERRQ(ierr);

    ierr = VecDestroy(&x);CHKERRQ(ierr);
    ierr = VecDestroy(&r);CHKERRQ(ierr);
    ierr = DMDestroy(&da);CHKERRQ(ierr);

    ierr = PetscFinalize();
    PetscFunctionReturn(0);
}

它可以通过以下方式编译:

mpicc main3.c -o main3 -O3 -I/.../petsc-3.6.4/include -I/.../petsc-3.6.4/linux-gnu-c-debug/include -L/.../petsc-3.6.4/linux-gnu-c-debug/lib -lpetsc

然后运行:

mpirun -np 4 main3 -da_grid_x 100 -da_grid_y 80