Значения vec PETSC DMDA назначаются на место проведения
Недавно я начал изучать 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?
Решение
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, в разделе 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);
Другие решения
Других решений пока нет …