c++ - 将PETSC DMDA vec值分配给调整位置

我最近开始学习PETSc，在尝试完成一些简单任务时遇到了问题。此代码有什么问题:

static char help[] = "Test  2d DMDAs Vecs.\n\n";
#include <petscdm.h>
#include <petscdmda.h>
#include <petscsys.h>
PetscReal process_value(int rank, int i) {
    return i*PetscPowScalar(10,rank*2);
}
int main(int argc,char **argv) {
  PetscErrorCode   ierr;
  PetscMPIInt      rank;
  PetscInt         M = -5,N = -3;
  DM               da;
  Vec              local,global;
  ierr = PetscInitialize(&argc,&argv,(char*)0,help);CHKERRQ(ierr);
  ierr  = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
  ierr = DMDACreate2d(PETSC_COMM_WORLD , DM_BOUNDARY_NONE , DM_BOUNDARY_NONE
    , DMDA_STENCIL_BOX , M , N , PETSC_DECIDE, PETSC_DECIDE
    ,   1 , 1 , NULL ,  NULL ,  &da); CHKERRQ(ierr);
  ierr = DMCreateGlobalVector(da,&global);CHKERRQ(ierr);
  ierr = DMCreateLocalVector(da,&local);CHKERRQ(ierr);
  {
      PetscInt       v,i, j,   xm, ym, xs, ys;
      PetscScalar    **array;
      ierr = DMDAGetCorners(da, &xs, &ys, 0, &xm, &ym, 0); CHKERRQ(ierr);
      PetscSynchronizedPrintf(PETSC_COMM_WORLD,"%d:xs=%d\txm=%d\tys=%d\tym=%d\n",rank,xs,xm,ys,ym);
      PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
      ierr = DMDAVecGetArray(da, global, &array); CHKERRQ(ierr);
    v=0;
      for (j = ys; j < ys + ym; j++) {
      for (i = xs; i < xs + xm; i++) {
          array[j][i] = process_value(rank,v+=1);
      }
      }
    ierr = DMDAVecRestoreArray(da, global, &array); CHKERRQ(ierr);
  }
    ierr = VecView(global,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  ierr = VecDestroy(&local);CHKERRQ(ierr);
  ierr = VecDestroy(&global);CHKERRQ(ierr);
  ierr = DMDestroy(&da);CHKERRQ(ierr);
  ierr = PetscFinalize();
  return 0;
}

它用过程等级标记的数量填充小阵列。
成功编译并链接后，将得到以下结果:

> mpiexec -n 2 ./problem
0:xs=0  xm=3    ys=0    ym=3
1:xs=3  xm=2    ys=0    ym=3
Vec Object: 2 MPI processes
  type: mpi
Vec Object:Vec_0x84000004_0 2 MPI processes
  type: mpi
Process [0]
1.
2.
3.
100.
200.
4.
5.
6.
300.
Process [1]
400.
7.
8.
9.
500.
600.
>

VecView显示流程已写入属于其他流程的位置。哪里有错误？ DMDAVecGetArray / DMDAVecRestoreArray给出错误的数组，或者VecView不适合查看从DM对象获得的Vec？

最佳答案

DMDAVecGetArray()和DMDAVecRestoreArray()可以正常工作。实际上，您正在处理几年前mailing list of petsc中描述的VecView()功能。

实际上，VecView()使用自然顺序从DMDA打印 vector 。

proc0   proc1
1  2  | 3  4
5  6  | 7  8
___________
9  10 | 11 12
13 14 | 15 16
proc2   proc3

在有关结构化网格的第2.5节documentation of Petsc中，强调了自然排序和Petsc排序之间的区别。 Petsc的订购如下所示:

proc0   proc1
1  2  | 5  6
3  4  | 7  8
___________
9  10 | 13 14
11 12 | 15 16
proc2   proc3

如线程VecView doesn't work properly with DA global vectors，中所示，仍然可以使用Petsc的顺序通过以下方式打印 vector :

PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_NATIVE);
VecView(global,PETSC_VIEWER_STDOUT_WORLD);
PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);

让我们仔细看看Petsc的来源，看看它是如何工作的。 DMDA vector 的VecView()操作被重载，因为调用了DMCreateGlobalVector_DA()函数(请参阅dadist.c):新方法是gr2.c中的VecView_MPI_DA()。毫不奇怪，它调用一个函数DMDACreateNaturalVector()，然后使用本机VecView()打印自然 vector 。如果使用格式PETSC_VIEWER_NATIVE，则vector interface调用*vec->ops->viewnative操作，该操作可能指向pdvec.c中的本机VecView()函数VecView_MPI_ASCII()。这解释了VecView用于DMDA vector 的奇怪(但非常实用!)行为。

如果希望保持自然顺序，则可以使用以下命令来删除打印出的无意义的Process[0]...Process[3]:

PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_COMMON);

关于c++ - 将PETSC DMDA vec值分配给调整位置，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/38794889/