Dear Petsc-Team,
Since Petsc-3.18 our code no longer runs successfully with periodic
geometries. Since this version, the call to DMGetPeriodicity has changed
arguments, but even after adapting our code, for some reason calling
DMLocalizeCoordinates, at least the way we do, no longer works. To
illustrate, I have made the following example, please find it attached.
The example creates a box mesh (10 x 10 x 10 cells) with unit size,
partitions this and distributes this. Then, the co-ordinates are
localized. The vertices of cell 9 (which on a 1 processor execution of
the code lies on a periodic domain) are printed before the partitioning
and after.
For PETSc versions <3.18, the output is:
Before partitioning/distributing/localizing:
Cell 9 has Vecclosuresize 24 and vertices:
9.000000e-01 0.000000e+00 0.000000e+00
9.000000e-01 1.000000e-01 0.000000e+00
0.000000e+00 1.000000e-01 0.000000e+00
0.000000e+00 0.000000e+00 0.000000e+00
9.000000e-01 0.000000e+00 1.000000e-01
0.000000e+00 0.000000e+00 1.000000e-01
0.000000e+00 1.000000e-01 1.000000e-01
9.000000e-01 1.000000e-01 1.000000e-01
After partitioning/distributing/localizing:
Cell 9 has Vecclosuresize 48 and vertices:
9.000000e-01 0.000000e+00 0.000000e+00
9.000000e-01 1.000000e-01 0.000000e+00
1.000000e+00 1.000000e-01 0.000000e+00
1.000000e+00 0.000000e+00 0.000000e+00
9.000000e-01 0.000000e+00 1.000000e-01
1.000000e+00 0.000000e+00 1.000000e-01
1.000000e+00 1.000000e-01 1.000000e-01
9.000000e-01 1.000000e-01 1.000000e-01
(the 'far away' X vertices are correctly put to 1.0)
But, with PETSc versions >= 3.18, we get the following output:
Before partitioning/distributing/localizing:
Cell 9 has Vecclosuresize 24 and vertices:
9.000000e-01 0.000000e+00 0.000000e+00
9.000000e-01 1.000000e-01 0.000000e+00
0.000000e+00 1.000000e-01 0.000000e+00
0.000000e+00 0.000000e+00 0.000000e+00
9.000000e-01 0.000000e+00 1.000000e-01
0.000000e+00 0.000000e+00 1.000000e-01
0.000000e+00 1.000000e-01 1.000000e-01
9.000000e-01 1.000000e-01 1.000000e-01
After partitioning/distributing/localizing:
Cell 9 has Vecclosuresize 24 and vertices:
9.000000e-01 0.000000e+00 0.000000e+00
9.000000e-01 1.000000e-01 0.000000e+00
0.000000e+00 1.000000e-01 0.000000e+00
0.000000e+00 0.000000e+00 0.000000e+00
9.000000e-01 0.000000e+00 1.000000e-01
0.000000e+00 0.000000e+00 1.000000e-01
0.000000e+00 1.000000e-01 1.000000e-01
9.000000e-01 1.000000e-01 1.000000e-01
and both the vertices as well as the Vecclosuresize seem to be different
from the earlier PETSc versions. Is this a bug? Or are we calling
functions in the incorrect order?
Thanks, best regards, Berend.
#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>
PetscErrorCode PrintVerticesDM(PetscInt rank, DM dm)
{
PetscSection section;
Vec coordinates;
PetscScalar *VecClosure = NULL;
PetscInt VecClosureSize;
PetscInt LocalCellID = 9;
PetscErrorCode ierr;
PetscFunctionBegin;
if (rank > 0) PetscFunctionReturn(0);
ierr = DMGetCoordinatesLocal(dm, &coordinates);
CHKERRQ(ierr);
ierr = DMGetCoordinateSection(dm, §ion);
CHKERRQ(ierr);
ierr = DMPlexVecGetClosure(dm, section, coordinates, LocalCellID, &VecClosureSize, &VecClosure);
CHKERRQ(ierr);
printf("Cell %i has Vecclosuresize %i and vertices:\n", LocalCellID, VecClosureSize);
for (PetscInt v = 0; v < 8; v++)
{
printf("%e %e %e\n", VecClosure[3 * v + 0], VecClosure[3 * v + 1], VecClosure[3 * v + 2]);
}
ierr = DMPlexVecRestoreClosure(dm, section, coordinates, LocalCellID, &VecClosureSize, &VecClosure);
CHKERRQ(ierr);
PetscFunctionReturn(0);
}
int main(int argc, char **argv)
{
const PetscBool cellSimplex = PETSC_FALSE;
const PetscBool interpBefore = PETSC_TRUE;
const PetscInt overlap = 1;
const PetscInt dim = 3;
PetscInt iSize[3];
PetscScalar Xmin[3], Xmax[3];
DMBoundaryType MeshPeriodicity[3]; /* Periodicity of mesh in each direction */
DM dm = NULL, 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);
/* Print co-ordinates first time, in PETSC < 3.18, these are not yet localized */
printf("Before partitioning/distributing/localizing:\n");
ierr = PrintVerticesDM(rank, 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);
printf("After partitioning/distributing/localizing:\n");
ierr = PrintVerticesDM(rank, dm);
CHKERRQ(ierr);
if (dm != NULL)
{
ierr = DMDestroy(&dm);
CHKERRQ(ierr);
}
/* Clean-up */
PetscFinalize();
return (0);
}
