c++ - PETSC DMDA vec 值分配给 wrang 位置

Question

我最近开始学习 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？

Answer 1

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

documentation of Petsc 中强调了自然排序和 Petsc 排序之间的区别。 ，在关于结构化网格的 2.5 节中。以下是 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 的源代码，看看它是如何工作的。当调用 DMCreateGlobalVector_DA() 函数时，DMDA vector 的 VecView() 操作被重载(参见 dadist.c ):新方法是 VecView_MPI_DA()在gr2.c .不出所料，它调用函数 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);