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存储为向量？你可能对Petsc的分布式数组感兴趣，参见DMDACreate2d() using 2 dofs u和v。然后，函数DMDAVecGetArrayDOF()将向量从DMDA转换为数组。最后，您的代码行可以编写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