2

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

4

1 回答 1

1

DMDAVecGetArray()并且DMDAVecRestoreArray()工作得很好。实际上，您正在处理几年前在 petsc 的邮件列表中VecView()描述的功能。

照原样VecView()，使用自然排序从 DMDA 打印向量。 

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

自然排序和 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 不能与 DA 全局向量正常工作一样，仍然可以使用 Petsc 的排序来打印向量，方法是：

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

让我们仔细看看 Petsc 的源代码，看看它是如何工作的。VecView()DMDA 向量的操作在DMCreateGlobalVector_DA()函数被调用时被重载（参见dadist.c）：新方法VecView_MPI_DA()在gr2.c中。不出所料，它调用了一个函数DMDACreateNaturalVector()，然后使用 native 打印自然向量VecView()。如果使用格式PETSC_VIEWER_NATIVE，向量接口调用操作*vec->ops->viewnative，它可能指向pdvec.c中的本机VecView()函数。这解释了 VecView 对 DMDA 向量的奇怪（但非常实用！）行为。VecView_MPI_ASCII()

如果您希望保持自然顺序，您可以使用以下方法取消打印的无意义Process[0]...Process[3]：

PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_COMMON);
于 2016-09-06T20:37:31.123 回答