I'm trying to use PETSc's DMDA Vectors with 2 degrees of freedom and access them using struct, like in the manual.
However, when I try to use DMDAVecGetArray
even with one degree of freedom (like in example below) I get memory double free or corruption
error. When I replace DMDAVecGetArray
with VecGetArray
everything works just fine.
What causes this error?
Compile MWE with
> gcc -Wall -Wextra -O0 -g $(pkg-config --cflags petsc mpi) -o mwe mwe.c $(pkg-config --libs petsc mpi)
and run with
> mpiexec --host localhost:$(nproc) -n $(nproc) mwe
MWE:
#include <petscdm.h>
#include <petscdmda.h>
#include <petscvec.h>
int
main(int argc, char **argv)
{
PetscErrorCode ierr;
#define E(x) \
do { \
ierr = x; \
CHKERRQ(ierr); \
} while (0)
PetscInt N = 10;
Vec fg, f;
ierr = PetscInitialize(&argc, &argv, NULL, NULL);
if (ierr) {
return ierr;
}
DM domain;
E(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_GHOSTED, N, 1, 1, NULL, &domain));
E(DMSetUp(domain));
E(DMCreateGlobalVector(domain, &fg));
E(DMCreateLocalVector(domain, &f));
E(VecZeroEntries(fg));
E(VecAssemblyBegin(fg));
E(VecAssemblyBegin(f));
E(VecAssemblyEnd(fg));
E(VecAssemblyEnd(f));
PetscInt size;
E(VecGetSize(f, &size));
E(DMGlobalToLocal(domain, fg, INSERT_VALUES, f));
PetscScalar *farr;
E(VecGetArray(f, &farr)); // replace this
// E(DMDAVecGetArrayWrite(domain, f, &farr)); // with this
for (PetscInt i = 0; i < size; i++) {
farr[i] += 1;
}
E(VecRestoreArray(f, &farr)); // replace this
// E(DMDAVecRestoreArrayWrite(domain, f, &farr)); // with this
E(DMLocalToGlobal(domain, f, INSERT_VALUES, fg));
E(VecDestroy(&f));
E(VecDestroy(&fg));
E(DMDestroy(&domain));
E(PetscFinalize());
return ierr;
}
Debian bullseye; petsc-dev 3.14.4+dfsg1-1; libopenmpi-dev 4.1.0-7