Dear PETSc Team,
We are currently having a problem with the function DMSetGlobalSection -
this worked with v3.16.2 and 3.16.1, but doesn't seem to work anymore in
PETSc versions after 3.16.2.
I've attached an example which replicates the issue. It creates a box
with a random field is generated, saved and loaded. The code worked for
PETSc v3.16.1 and v3.16.2 (main branch) but fails for the last main (and
release) version v3.16.4.
The issue arises in DMSetGlobalSection(sdm,s) (line 160) as it detects
the cloned DM sdm as Null argument. I've also included the error message
at the bottom.
Are we doing something wrong? Or has the 'Cloning' of the dm changed
recently?
Thanks and best regards,
Berend
********************** Error message ***********************************
[0]PETSC ERROR: --------------------- Error Message
--------------------------------------------------------------
[0]PETSC ERROR: Null argument, when expecting valid pointer
[0]PETSC ERROR: Null Object: Parameter # 1
[0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.
[0]PETSC ERROR: Petsc Development GIT revision: v3.16.4-1031-g8387aa1ef1
GIT Date: 2022-02-21 23:19:21 -0600
[0]PETSC ERROR: ./bin/restart_periodic on a arch-linux-c-opt named
zeldovich by serbenlo Tue Feb 22 14:32:56 2022
[0]PETSC ERROR: Configure options --with-debugging=yes
--with-errorchecking=yes --with-clean --download-metis=yes
--download-parmetis=yes --download-hdf5 --download-p4est
--download-triangle --download-tetgen
--with-zlib-lib=/usr/lib/x86_64-linux-gnu/libz.a
--with-zlib-include=/usr/include --with-mpi=yes --with-mpi-dir=/usr
--with-mpiexec=/usr/bin/mpiexec
[0]PETSC ERROR: #1 PetscSectionGetDof() at
/usr/local/petsc_main/src/vec/is/section/interface/section.c:807
[0]PETSC ERROR: #2 DMDefaultSectionCheckConsistency_Internal() at
/usr/local/petsc_main/src/dm/interface/dm.c:4520
[0]PETSC ERROR: #3 DMSetGlobalSection() at
/usr/local/petsc_main/src/dm/interface/dm.c:4614
[0]PETSC ERROR: #4 main() at
/home/serbenlo/DATA/Arbeit/OvGU_PostDoc_2021/Projects/MF_Restart/Periodic_DM/restart-periodic/src/main.c:161
[0]PETSC ERROR: No PETSc Option Table entries
[0]PETSC ERROR: ----------------End of Error Message -------send entire
error message to petsc-ma...@mcs.anl.gov----------
application called MPI_Abort(MPI_COMM_WORLD, 85) - process 0
#include "petsc.h"
#include "petscdmplex.h"
#include "petscdm.h"
#include "petscdmlabel.h"
#include "petscviewerhdf5.h"
#include "petscds.h"
#include "petscvec.h"
#include <time.h>
#include <stdlib.h>
int main(int argc, char **argv)
{
/**** Generate Periodic DM ******/
const PetscBool cellSimplex = PETSC_FALSE;
const PetscBool interpBefore = PETSC_TRUE;
PetscInt dim = 3;
PetscInt iSize[3];
PetscInt i, numPoints;
PetscInt overlap = 1, numFields = 1;
PetscScalar Xmin[3], Xmax[3];
PetscScalar *xVecArray;
DMBoundaryType MeshPeriodicity[3]; /* Periodicity of mesh in each direction */
DM dm, pdm = NULL;
PetscPartitioner part;
PetscErrorCode ierr;
srand(time(NULL)); /* Random values initialization */
ierr = PetscInitialize(&argc, &argv, NULL, NULL);
CHKERRQ(ierr);
// Number of cells in each direction
iSize[0] = 10;
iSize[1] = 10;
iSize[2] = 10;
// Dimensions
Xmin[0] = 0.0;
Xmin[1] = 0.0;
Xmin[2] = 0.0;
Xmax[0] = 1.0;
Xmax[1] = 1.0;
Xmax[2] = 1.0;
// Periodicity
MeshPeriodicity[0] = DM_BOUNDARY_NONE;
MeshPeriodicity[1] = DM_BOUNDARY_NONE;
MeshPeriodicity[2] = DM_BOUNDARY_NONE;
// Create Mesh
ierr = DMPlexCreateBoxMesh(PETSC_COMM_WORLD,dim,cellSimplex,iSize,Xmin,Xmax,MeshPeriodicity,interpBefore,&dm);
CHKERRQ(ierr);
// Create Partition
ierr = DMPlexGetPartitioner(dm, &part);
CHKERRQ(ierr);
if (part)
{
ierr = PetscPartitionerSetType(part, PETSCPARTITIONERPARMETIS);
CHKERRQ(ierr);
ierr = PetscPartitionerSetFromOptions(part);
CHKERRQ(ierr);
}
// Distribute
ierr = DMPlexDistribute(dm, overlap, NULL, &pdm);
CHKERRQ(ierr);
if (pdm)
{
ierr = DMDestroy(&dm);
CHKERRQ(ierr);
dm = pdm;
ierr = PetscObjectSetName((PetscObject) dm, "Distributed Mesh");
CHKERRQ(ierr);
}
ierr = DMLocalizeCoordinates(dm);
CHKERRQ(ierr);
ierr = DMSetFromOptions(dm);
CHKERRQ(ierr);
// Create Section on DM
PetscSection section;
PetscInt *numComp;
PetscInt *numDof;
ierr = PetscCalloc2(numFields, &numComp, numFields * (dim + 1), &numDof);
CHKERRQ(ierr);
for (i = 0; i < numFields; i++)
{
numComp[i] = 1;
numDof[i * (dim + 1) + dim] = 1;
}
ierr = DMSetNumFields(dm, numFields);
CHKERRQ(ierr);
ierr = DMPlexCreateSection(dm, NULL, numComp, numDof, 0, NULL, NULL, NULL, NULL, §ion);
CHKERRQ(ierr);
// Push section onto the DM
ierr = DMSetLocalSection(dm, section);
CHKERRQ(ierr);
ierr = PetscFree2(numComp, numDof);
CHKERRQ(ierr);
ierr = PetscSectionDestroy(§ion);
CHKERRQ(ierr);
/**** Create Vector Field => xGlobalVector *********/
Vec xGlobalVector;
ierr = DMCreateGlobalVector(dm, &xGlobalVector);
CHKERRQ(ierr);
// Initialize the Vec to 0
ierr = VecSet(xGlobalVector,0.0);
CHKERRQ(ierr);
// Fill the Vector with random values
ierr = VecGetArray(xGlobalVector, &xVecArray);
CHKERRQ(ierr);
ierr = VecGetLocalSize(xGlobalVector, &numPoints);
CHKERRQ(ierr);
for (i = 0; i < numPoints; i++) /* Loop over all internal cells */
{
xVecArray[i] = rand() % 20;
}
ierr = VecRestoreArray(xGlobalVector, &xVecArray);
CHKERRQ(ierr);
/**** Write DM + associated Vec ****/
PetscViewer H5Viewer;
ierr = PetscViewerHDF5Open(PETSC_COMM_WORLD, "periodic_box.h5", FILE_MODE_WRITE, &H5Viewer);
CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject)dm, "plexA");
CHKERRQ(ierr);
ierr = PetscViewerPushFormat(H5Viewer, PETSC_VIEWER_HDF5_PETSC);
CHKERRQ(ierr);
ierr = DMPlexTopologyView(dm, H5Viewer);
CHKERRQ(ierr);
ierr = DMPlexLabelsView(dm, H5Viewer);
CHKERRQ(ierr);
ierr = DMPlexCoordinatesView(dm, H5Viewer);
CHKERRQ(ierr);
ierr = PetscViewerPopFormat(H5Viewer);
CHKERRQ(ierr);
DM sdm;
PetscSection s;
ierr = DMClone(dm, &sdm);
CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject)sdm, "dmA");
CHKERRQ(ierr);
ierr = DMGetGlobalSection(dm, &s);
CHKERRQ(ierr);
ierr = DMSetGlobalSection(sdm, s);
CHKERRQ(ierr);
ierr = DMPlexSectionView(dm, H5Viewer, sdm);
CHKERRQ(ierr);
Vec vec;
PetscScalar *array;
ierr = DMGetGlobalVector(sdm, &vec);
CHKERRQ(ierr);
/*** Fill the vector xGlobalVector ***/
ierr = VecGetArray(vec, &array);
CHKERRQ(ierr);
ierr = VecGetLocalSize(xGlobalVector, &numPoints);
CHKERRQ(ierr);
ierr = VecGetArray(xGlobalVector, &xVecArray);
CHKERRQ(ierr);
for (i = 0; i < numPoints; i++) /* Loop over all internal cells */
{
array[i] = xVecArray[i];
}
ierr = VecRestoreArray(vec, &array);
CHKERRQ(ierr);
ierr = VecRestoreArray(xGlobalVector, &xVecArray);
CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject)vec, "vecA");
CHKERRQ(ierr);
ierr = DMPlexGlobalVectorView(dm, H5Viewer, sdm, vec);
CHKERRQ(ierr);
ierr = PetscViewerDestroy(&H5Viewer);
CHKERRQ(ierr);
/*** End of writing ****/
/*** Wait for all CPUs to end ***/
ierr = MPI_Barrier(PETSC_COMM_WORLD);
CHKERRQ(ierr);
/*** Load ***/
ierr = DMDestroy(&dm);
CHKERRQ(ierr);
ierr = DMDestroy(&sdm);
CHKERRQ(ierr);
Vec Restart_xGlobalVector;
PetscSF sfO, sf, sfDist, globalDataSF, localDataSF;
DM distributedDM;
ierr = PetscViewerHDF5Open(PETSC_COMM_WORLD, "periodic_box.h5", FILE_MODE_READ, &H5Viewer);
CHKERRQ(ierr);
ierr = DMCreate(PETSC_COMM_WORLD, &dm);
CHKERRQ(ierr);
ierr = DMSetType(dm, DMPLEX);
CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject)dm, "plexA");
CHKERRQ(ierr);
ierr = PetscViewerPushFormat(H5Viewer, PETSC_VIEWER_HDF5_PETSC);
CHKERRQ(ierr);
ierr = DMPlexTopologyLoad(dm, H5Viewer, &sfO);
CHKERRQ(ierr);
ierr = DMPlexLabelsLoad(dm, H5Viewer, sfO);
CHKERRQ(ierr);
ierr = DMPlexCoordinatesLoad(dm, H5Viewer, sfO);
CHKERRQ(ierr);
ierr = PetscViewerPopFormat(H5Viewer);
CHKERRQ(ierr);
ierr = DMPlexDistribute(dm, overlap, &sfDist, &distributedDM);
CHKERRQ(ierr);
if (distributedDM) {
ierr = DMDestroy(&dm);
CHKERRQ(ierr);
dm = distributedDM;
ierr = PetscObjectSetName((PetscObject)dm, "plexA");
CHKERRQ(ierr);
}
ierr = PetscSFCompose(sfO, sfDist, &sf);
CHKERRQ(ierr);
ierr = PetscSFDestroy(&sfO);
CHKERRQ(ierr);
ierr = PetscSFDestroy(&sfDist);
CHKERRQ(ierr);
ierr = DMClone(dm, &sdm);
CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject)sdm, "dmA");
CHKERRQ(ierr);
ierr = DMPlexSectionLoad(dm, H5Viewer, sdm, sf, &globalDataSF, &localDataSF);
CHKERRQ(ierr);
/*** Load the Vector to Restart_xGlobalVector ***/
ierr = DMGetGlobalVector(sdm, &Restart_xGlobalVector);
CHKERRQ(ierr);
ierr = VecSet(Restart_xGlobalVector,0.0);
CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject)Restart_xGlobalVector, "vecA");
CHKERRQ(ierr);
ierr = DMPlexGlobalVectorLoad(dm, H5Viewer, sdm,globalDataSF,Restart_xGlobalVector);
CHKERRQ(ierr);
ierr = PetscViewerDestroy(&H5Viewer);
CHKERRQ(ierr);
PetscEnd();
}
