Dear PETSc team and users,
I have asked a few times about this before, but we haven't really gotten
this to work yet.
In our code, we use the DMPlex framework and are also interested in
periodic geometries.
As our simulations typically require many time-steps, we would like to
be able to save the DM to file and to read it again to resume the
simulation (a restart).
Although this works for a non-periodic DM, we haven't been able to get
this to work for a periodic one. To illustrate this, I have made a
working example, consisting of 2 files, createandwrite.c and
readandcreate.c. I have attached these 2 working examples. We are using
Petsc-3.18.2.
In the first file (createandwrite.c) a DMPlex is created and written to
a file. Periodicity is activated on lines 52-55 of the code.
In the second file (readandcreate.c) a DMPlex is read from the file.
When a periodic DM is read, this does not work. Also, trying to
'enforce' periodicity, lines 55 - 66, does not work if the number of
processes is larger than 1 - the code "hangs" without producing an error.
Could you indicate what I am missing? I have really tried many different
options, without finding a solution.
Many thanks and kind regards,
Berend.
// Create a DM like MF3, and write it to disc with HDF5.
#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>
#include <assert.h>
int main(int argc, char **argv)
{
const PetscBool cellSimplex = PETSC_FALSE;
const PetscBool interpBefore = PETSC_TRUE;
const PetscInt overlap = 1, numFields = 4; /* Overlap and number of fields */
PetscInt dim = 3;
PetscInt iSize[3];
PetscInt i, numPoints;
PetscScalar Xmin[3], Xmax[3];
PetscScalar *xVecArray;
DMBoundaryType MeshPeriodicity[3]; /* Periodicity of mesh in each direction */
DM dm, pdm = NULL;
PetscPartitioner partitioner;
PetscErrorCode ierr;
PetscInt rank, psize;
ierr = PetscInitialize(&argc, &argv, NULL, NULL);
CHKERRQ(ierr);
ierr = MPI_Comm_rank(MPI_COMM_WORLD, &rank);
CHKERRQ(ierr);
ierr = MPI_Comm_size(MPI_COMM_WORLD, &psize);
CHKERRQ(ierr);
// Number of cells in each direction
iSize[0] = 10;
iSize[1] = 10;
iSize[2] = 10;
// Physical 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_PERIODIC;
MeshPeriodicity[1] = DM_BOUNDARY_PERIODIC;
MeshPeriodicity[2] = DM_BOUNDARY_PERIODIC;
// MeshPeriodicity[0] = DM_BOUNDARY_NONE;
// MeshPeriodicity[1] = DM_BOUNDARY_NONE;
// MeshPeriodicity[2] = DM_BOUNDARY_NONE;
// Create Mesh with DMPlex
ierr = DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, cellSimplex, iSize, Xmin, Xmax, MeshPeriodicity, interpBefore, &dm);
CHKERRQ(ierr);
// Partition the mesh with Parmetis
ierr = DMPlexGetPartitioner(dm, &partitioner);
CHKERRQ(ierr);
if (partitioner)
{
ierr = PetscPartitionerSetType(partitioner, PETSCPARTITIONERPARMETIS);
CHKERRQ(ierr);
ierr = PetscPartitionerSetFromOptions(partitioner);
CHKERRQ(ierr);
}
ierr = DMSetBasicAdjacency(dm, PETSC_TRUE, PETSC_TRUE);
CHKERRQ(ierr);
// Distribute the mesh
ierr = DMPlexDistribute(dm, overlap, NULL, &pdm);
CHKERRQ(ierr);
if (pdm)
{
ierr = DMDestroy(&dm);
CHKERRQ(ierr);
dm = pdm;
}
ierr = PetscObjectSetName((PetscObject) dm, "DM");
CHKERRQ(ierr);
ierr = DMPlexDistributeSetDefault(dm, PETSC_FALSE);
CHKERRQ(ierr);
/* Localize the co-ordinates, required for periodic domains */
ierr = DMLocalizeCoordinates(dm);
CHKERRQ(ierr);
ierr = DMSetFromOptions(dm);
CHKERRQ(ierr);
// Create Section on DM with fields and degrees of Freedom
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);
/* Now free/release the variables for the section and destroy the section */
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 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] = (PetscScalar) i;
}
ierr = VecRestoreArray(xGlobalVector, &xVecArray);
CHKERRQ(ierr);
/*****************************************************************************/
/* We are now in the same situation as in MF3. */
/*****************************************************************************/
/**** Write DM + associated Vec ****/
PetscViewer H5Viewer;
ierr = PetscViewerHDF5Open(PETSC_COMM_WORLD, "mesh.h5", FILE_MODE_WRITE, &H5Viewer);
CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject) dm, "DM");
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 dmCopy;
PetscSection sectionCopy;
ierr = DMClone(dm, &dmCopy);
CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject) dmCopy, "DM");
CHKERRQ(ierr);
ierr = DMGetLocalSection(dm, §ionCopy);
CHKERRQ(ierr);
ierr = DMSetLocalSection(dmCopy, sectionCopy);
CHKERRQ(ierr);
/* Write the section to the file */
ierr = DMPlexSectionView(dm, H5Viewer, dmCopy);
CHKERRQ(ierr);
Vec vectorCopy;
PetscScalar *array;
PetscInt numPointsCopy;
ierr = DMGetGlobalVector(dmCopy, &vectorCopy);
CHKERRQ(ierr);
/*** We have to copy the vector into the new vector ... ***/
ierr = VecGetArray(vectorCopy, &array);
CHKERRQ(ierr);
ierr = VecGetLocalSize(xGlobalVector, &numPoints);
CHKERRQ(ierr);
ierr = VecGetLocalSize(vectorCopy, &numPointsCopy);
CHKERRQ(ierr);
assert(numPoints == numPointsCopy);
ierr = VecGetArray(xGlobalVector, &xVecArray);
CHKERRQ(ierr);
for (i = 0; i < numPoints; i++) /* Loop over all internal cells */
{
array[i] = xVecArray[i];
}
ierr = VecRestoreArray(vectorCopy, &array);
CHKERRQ(ierr);
ierr = VecRestoreArray(xGlobalVector, &xVecArray);
CHKERRQ(ierr);
// Output information on mesh which has just been written:
PetscInt MaxConeSize, MaxSupportSize, CellHeight;
const PetscScalar *MaxCell, *Lstart, *L;
PetscInt cStart, cEnd;
ierr = DMPlexGetMaxSizes(dmCopy, &MaxConeSize, &MaxSupportSize);
CHKERRQ(ierr);
ierr = DMPlexGetVTKCellHeight(dmCopy, &CellHeight);
CHKERRQ(ierr);
ierr = DMPlexGetHeightStratum(dmCopy, CellHeight, &cStart, &cEnd);
CHKERRQ(ierr);
ierr = DMGetDimension(dm, &dim);
CHKERRQ(ierr);
ierr = DMGetPeriodicity(dmCopy, &MaxCell, &Lstart, &L);
CHKERRQ(ierr);
for (i = 0; i < psize; i++)
{
if (i == rank)
{
PetscPrintf(PETSC_COMM_SELF, "Processor %i from %i; DM Dimension %i, max Cone Size %i, MaxSupportSize %i, CellHeight %i, cStart: %i -> Cend %i ", rank,
psize, dim, MaxConeSize, MaxSupportSize, CellHeight, cStart, cEnd);
if (L != NULL)
{
PetscPrintf(PETSC_COMM_SELF, " PeriodicLengths: ");
if (L[0])
{
PetscPrintf(PETSC_COMM_SELF, " X: %e,%e,%e", MaxCell[0], Lstart[0], L[0]);
}
if (L[1])
{
PetscPrintf(PETSC_COMM_SELF, " Y: %e,%e,%e", MaxCell[1], Lstart[1], L[1]);
}
if (L[2])
{
PetscPrintf(PETSC_COMM_SELF, " Z: %e,%e,%e", MaxCell[2], Lstart[2], L[2]);
}
PetscPrintf(PETSC_COMM_SELF, "\n ");
}
else
{
PetscPrintf(PETSC_COMM_SELF, " Mesh not periodic \n");
}
}
ierr = MPI_Barrier(PETSC_COMM_WORLD);
CHKERRQ(ierr);
}
ierr = PetscObjectSetName((PetscObject) vectorCopy, "DataVector");
CHKERRQ(ierr);
/* Write the vector to the file */
ierr = DMPlexGlobalVectorView(dm, H5Viewer, dmCopy, vectorCopy);
CHKERRQ(ierr);
/* Close the file */
ierr = PetscViewerDestroy(&H5Viewer);
CHKERRQ(ierr);
/*** End of writing ****/
PetscFinalize();
return (0);
}
// Read an existing DM and vector
#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>
#include <assert.h>
int main(int argc, char **argv)
{
PetscInt i;
const PetscInt overlap = 1;
DM dm, dmCopy, distributedDM;
PetscPartitioner partitioner;
PetscErrorCode ierr;
Vec Restart_xGlobalVector;
PetscSF sfO, sf, sfDist, globalDataSF, localDataSF;
PetscViewer H5Viewer;
PetscInt rank, psize;
ierr = PetscInitialize(&argc, &argv, NULL, NULL);
CHKERRQ(ierr);
ierr = MPI_Comm_rank(MPI_COMM_WORLD, &rank);
CHKERRQ(ierr);
ierr = MPI_Comm_size(MPI_COMM_WORLD, &psize);
CHKERRQ(ierr);
/* Open the file */
ierr = PetscViewerHDF5Open(PETSC_COMM_WORLD, "mesh.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, "DM");
CHKERRQ(ierr);
ierr = PetscViewerPushFormat(H5Viewer, PETSC_VIEWER_HDF5_PETSC);
CHKERRQ(ierr);
/* Load the mesh in the same order from the HDF5 file */
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);
// /* Try to FORCE periodicity: */
// {
// const PetscScalar maxCell[] = {0.11, 0.11, 0.11};
// const PetscScalar Lstart[] = {0.0, 0.0, 0.0};
// const PetscScalar L[] = {1.0, 1.0, 1.0};
// ierr = DMSetPeriodicity(dm, maxCell, Lstart, L);
// CHKERRQ(ierr);
// }
// Partition the mesh with Parmetis
ierr = DMPlexGetPartitioner(dm, &partitioner);
CHKERRQ(ierr);
if (partitioner)
{
ierr = PetscPartitionerSetType(partitioner, PETSCPARTITIONERPARMETIS);
CHKERRQ(ierr);
ierr = PetscPartitionerSetFromOptions(partitioner);
CHKERRQ(ierr);
}
/* Localize the co-ordinates, required for periodic domains */
ierr = DMLocalizeCoordinates(dm);
CHKERRQ(ierr);
ierr = DMSetFromOptions(dm);
CHKERRQ(ierr);
ierr = DMPlexDistribute(dm, overlap, &sfDist, &distributedDM);
CHKERRQ(ierr);
if (distributedDM)
{
ierr = DMDestroy(&dm);
CHKERRQ(ierr);
dm = distributedDM;
ierr = PetscObjectSetName((PetscObject) dm, "DM");
CHKERRQ(ierr);
}
if (sfDist != NULL)
{
ierr = PetscSFCompose(sfO, sfDist, &sf);
CHKERRQ(ierr);
ierr = PetscSFDestroy(&sfO);
CHKERRQ(ierr);
ierr = PetscSFDestroy(&sfDist);
CHKERRQ(ierr);
}
else
{
sf = sfO;
}
ierr = DMPlexDistributeSetDefault(dm, PETSC_FALSE);
CHKERRQ(ierr);
ierr = DMClone(dm, &dmCopy);
CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject) dmCopy, "DM");
CHKERRQ(ierr);
ierr = DMPlexSectionLoad(dm, H5Viewer, dmCopy, sf, &globalDataSF, &localDataSF);
CHKERRQ(ierr);
/*** Load the Vector to Restart_xGlobalVector ***/
ierr = DMGetGlobalVector(dmCopy, &Restart_xGlobalVector);
CHKERRQ(ierr);
ierr = VecSet(Restart_xGlobalVector, 0.0);
CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject) Restart_xGlobalVector, "DataVector");
CHKERRQ(ierr);
ierr = DMPlexGlobalVectorLoad(dm, H5Viewer, dmCopy, globalDataSF, Restart_xGlobalVector);
CHKERRQ(ierr);
// Report on DM and vector read
PetscInt MaxConeSize, MaxSupportSize, CellHeight;
const PetscScalar *MaxCell, *Lstart, *L;
PetscInt cStart, cEnd;
PetscInt dim;
ierr = DMPlexGetMaxSizes(dmCopy, &MaxConeSize, &MaxSupportSize);
CHKERRQ(ierr);
ierr = DMPlexGetVTKCellHeight(dmCopy, &CellHeight);
CHKERRQ(ierr);
ierr = DMPlexGetHeightStratum(dmCopy, CellHeight, &cStart, &cEnd);
CHKERRQ(ierr);
ierr = DMGetDimension(dm, &dim);
CHKERRQ(ierr);
ierr = DMGetPeriodicity(dmCopy, &MaxCell, &Lstart, &L);
CHKERRQ(ierr);
// Output information on read mesh:
ierr = DMGetPeriodicity(dmCopy, &MaxCell, &Lstart, &L);
CHKERRQ(ierr);
for (i = 0; i < psize; i++)
{
if (i == rank)
{
PetscPrintf(PETSC_COMM_SELF, "Processor %i from %i; DM Dimension %i, max Cone Size %i, MaxSupportSize %i, CellHeight %i, cStart: %i -> Cend %i ", rank,
psize, dim, MaxConeSize, MaxSupportSize, CellHeight, cStart, cEnd);
if (L != NULL)
{
PetscPrintf(PETSC_COMM_SELF, " PeriodicLengths: ");
if (L[0])
{
PetscPrintf(PETSC_COMM_SELF, " X: %e,%e,%e", MaxCell[0], Lstart[0], L[0]);
}
if (L[1])
{
PetscPrintf(PETSC_COMM_SELF, " Y: %e,%e,%e", MaxCell[1], Lstart[1], L[1]);
}
if (L[2])
{
PetscPrintf(PETSC_COMM_SELF, " Z: %e,%e,%e", MaxCell[2], Lstart[2], L[2]);
}
PetscPrintf(PETSC_COMM_SELF, "\n ");
}
else
{
PetscPrintf(PETSC_COMM_SELF, " Mesh not periodic \n");
}
}
ierr = MPI_Barrier(PETSC_COMM_WORLD);
CHKERRQ(ierr);
}
// VecView(Restart_xGlobalVector, PETSC_VIEWER_STDOUT_WORLD);
/* Clean-up */
ierr = DMRestoreGlobalVector(dmCopy, &Restart_xGlobalVector);
CHKERRQ(ierr);
ierr = DMDestroy(&dm);
CHKERRQ(ierr);
ierr = DMDestroy(&dmCopy);
CHKERRQ(ierr);
ierr = PetscViewerDestroy(&H5Viewer);
CHKERRQ(ierr);
PetscFinalize();
return (0);
}
