Значения 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 не подходит для просмотра Vec, полученного из объекта DM?


c++ petsc
person avb1003    schedule 05.08.2016    source источник

Ответы (1)


arrow_upward
1
arrow_downward

DMDAVecGetArray() и DMDAVecRestoreArray() работают отлично. Действительно, вы имеете дело с функцией VecView(), описанной в список рассылки petsc несколько лет назад.

Как есть, 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(), а затем печатает естественный вектор, используя собственный VecView(). Если используется формат PETSC_VIEWER_NATIVE, векторный интерфейс вызывает операцию *vec->ops->viewnative, которая, вероятно, указывает на родную VecView() функцию VecView_MPI_ASCII() в pdvec.c. Это объясняет странное (но очень практичное!) поведение VecView для векторов DMDA.

Если вы хотите сохранить естественный порядок, вы можете удалить бессмысленные Process[0]...Process[3], которые печатаются, используя:

PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_COMMON);
person francis    schedule 06.09.2016