Hi Matthew,
I have joined the ex44.c example so you can see that we see.
The problem is about the "f" field number which I can't see any clue in
the documentation to what to put there... We are missing some
information maybe written elsewhere?
I have tried 3 different calls (lines 180-182) with comments aside:
//ierr = DMSetBasicAdjacency(ddm, PETSC_TRUE, PETSC_FALSE);
CHKERRQ(ierr); //This works: it give only 2 cells for the overlap
//ierr = DMSetAdjacency(ddm, 0, PETSC_TRUE, PETSC_FALSE); CHKERRQ(ierr);
//This works too, but why do I have to use field #0 ??? Is it a
convention? ierr = DMSetAdjacency(ddm, PETSC_DEFAULT, PETSC_TRUE,
PETSC_FALSE); CHKERRQ(ierr); //This does not works: it give 3 cells for
the overlap instead of 2.
Also, do you think your gitlab/knepley/fix-plex-g2n branch will find
it's way into next or main? ;)
Thanks,
Eric
On 2021-11-10 3:26 p.m., Matthew Knepley wrote:
Okay, so the PETSC_DEFAULT one is what is used for topology. There
are two separate uses for adjacency info.
First, for this kind of topological extension, and the second is for
calculating Jacobians. You should be able to see
the difference between FEM and FVM in just the case of a 2x2 mesh
where each process gets 1 cell. With FEM,
the overlap is the whole mesh, whereas with FVM it leaves out the
diagonal partner.
Thanks,
Matt
static char help[] = "The main goal of this code is to retrieve the original element numbers as found in the "
"initial partitions (sInitialPartition)... but after the call to DMPlexDistribute";
#include <petsc.h>
PetscReal sCoords2x5Mesh[18][2] = {
{0.00000000000000000e+00, 0.00000000000000000e+00},
{2.00000000000000000e+00, 0.00000000000000000e+00},
{0.00000000000000000e+00, 1.00000000000000000e+00},
{2.00000000000000000e+00, 1.00000000000000000e+00},
{9.99999999997387978e-01, 0.00000000000000000e+00},
{9.99999999997387978e-01, 1.00000000000000000e+00},
{0.00000000000000000e+00, 2.00000000000000011e-01},
{0.00000000000000000e+00, 4.00000000000000022e-01},
{0.00000000000000000e+00, 5.99999999999999978e-01},
{0.00000000000000000e+00, 8.00000000000000044e-01},
{2.00000000000000000e+00, 2.00000000000000011e-01},
{2.00000000000000000e+00, 4.00000000000000022e-01},
{2.00000000000000000e+00, 5.99999999999999978e-01},
{2.00000000000000000e+00, 8.00000000000000044e-01},
{9.99999999997387756e-01, 2.00000000000000011e-01},
{9.99999999997387978e-01, 4.00000000000000022e-01},
{9.99999999997387978e-01, 6.00000000000000089e-01},
{9.99999999997388089e-01, 8.00000000000000044e-01}};
//Connectivity of a 2x5 rectangular mesh of quads :
const PetscInt sConnectivity2x5Mesh[10][4] = {
{0,4,14,6},
{6,14,15,7},
{7,15,16,8},
{8,16,17,9},
{9,17,5,2},
{4,1,10,14},
{14,10,11,15},
{15,11,12,16},
{16,12,13,17},
{17,13,3,5}};
const PetscInt sInitialPartition2x5Mesh[2][5] = {
{0,2,4,6,8},
{1,3,5,7,9}
};
const PetscInt sNLoclCells2x5Mesh = 5;
const PetscInt sNGlobVerts2x5Mesh = 18;
int main(int argc, char **argv)
{
const PetscInt Nc = sNLoclCells2x5Mesh; //Same on each rank for this example...
const PetscInt Nv = sNGlobVerts2x5Mesh;
const PetscInt* InitPartForRank[2] = {&sInitialPartition2x5Mesh[0][0],
&sInitialPartition2x5Mesh[1][0]};
const PetscInt (*Conn)[4] = sConnectivity2x5Mesh;
const PetscInt Ncor = 4;
const PetscInt dim = 2;
DM dm, idm, ddm;
PetscSF sfVert, sfMig, sfPart;
PetscPartitioner part;
PetscSection s;
PetscInt *cells, c;
PetscMPIInt size, rank;
PetscBool box = PETSC_FALSE, field = PETSC_FALSE;
PetscErrorCode ierr;
ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr;
ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size);CHKERRMPI(ierr);
ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);CHKERRMPI(ierr);
if (size != 2) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "This is a 2 processors example only");
ierr = PetscOptionsGetBool(NULL, NULL, "-box", &box, NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetBool(NULL, NULL, "-field", &field, NULL);CHKERRQ(ierr);
ierr = DMPlexCreate(PETSC_COMM_WORLD, &dm);CHKERRQ(ierr);
if (box) {
ierr = DMSetType(dm, DMPLEX);CHKERRQ(ierr);
ierr = DMSetFromOptions(dm);CHKERRQ(ierr);
} else {
ierr = PetscMalloc1(Nc * Ncor, &cells);CHKERRQ(ierr);
for (c = 0; c < Nc; ++c) {
PetscInt cell = (InitPartForRank[rank])[c], cor;
for (cor = 0; cor < Ncor; ++cor) {
cells[c*Ncor + cor] = Conn[cell][cor];
}
}
ierr = DMSetDimension(dm, dim);CHKERRQ(ierr);
ierr = DMPlexBuildFromCellListParallel(dm, Nc, PETSC_DECIDE, Nv, Ncor, cells, &sfVert);CHKERRQ(ierr);
//ierr = DMPlexBuildCoordinatesFromCellListParallel(dm, dim, sfVert, coords);CHKERRQ(ierr);
ierr = PetscSFDestroy(&sfVert);CHKERRQ(ierr);
ierr = PetscFree(cells);CHKERRQ(ierr);
ierr = DMPlexInterpolate(dm, &idm);CHKERRQ(ierr);
ierr = DMDestroy(&dm);CHKERRQ(ierr);
dm = idm;
}
ierr = DMSetUseNatural(dm, PETSC_TRUE);CHKERRQ(ierr);
ierr = DMViewFromOptions(dm, NULL, "-dm_view");CHKERRQ(ierr);
if (field) {
const PetscInt Nf = 1;
const PetscInt numComp[1] = {1};
const PetscInt numDof[3] = {0, 0, 1};
const PetscInt numBC = 0;
ierr = DMSetNumFields(dm, Nf); CHKERRQ(ierr);
ierr = DMPlexCreateSection(dm, NULL, numComp, numDof, numBC, NULL, NULL, NULL, NULL, &s); CHKERRQ(ierr);
ierr = DMSetLocalSection(dm, s); CHKERRQ(ierr);
ierr = PetscSectionView(s, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
ierr = PetscSectionDestroy(&s);CHKERRQ(ierr);
}
ierr = DMPlexGetPartitioner(dm, &part);CHKERRQ(ierr);
ierr = PetscPartitionerSetFromOptions(part);CHKERRQ(ierr);
ierr = DMPlexDistribute(dm, 0, &sfMig, &ddm); CHKERRQ(ierr);
ierr = PetscSFView(sfMig, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
ierr = PetscSFCreateInverseSF(sfMig, &sfPart); CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject) sfPart, "Inverse Migration SF");CHKERRQ(ierr);
ierr = PetscSFView(sfPart, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
Vec lNatVec ;
Vec lGlobalVec;
ierr = DMGetGlobalVector(dm,&lNatVec); CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject) lNatVec, "Natural Vector (initial partition)");CHKERRQ(ierr);
PetscScalar *lNatVecArray;
ierr = VecGetArray(lNatVec, &lNatVecArray); CHKERRQ(ierr);
//Copying the initial partition into the "natural" vector:
for (c = 0; c < Nc; ++c) {
lNatVecArray[c] = (InitPartForRank[rank])[c];
}
ierr = VecRestoreArray(lNatVec, &lNatVecArray); CHKERRQ(ierr);
ierr = DMGetGlobalVector(ddm,&lGlobalVec); CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject) lGlobalVec, "Global Vector (reordered element numbers in the petsc distributed order)");CHKERRQ(ierr);
ierr = VecZeroEntries(lGlobalVec); CHKERRQ(ierr);
// The call to DMPlexNaturalToGlobalBegin/End does not produce our expected result...
// In lGlobalVec, we expect to have:
/*
* Process [0]
* 2.
* 4.
* 8.
* 3.
* 9.
* Process [1]
* 1.
* 5.
* 7.
* 0.
* 6.
*/
ierr = DMPlexNaturalToGlobalBegin(ddm, lNatVec, lGlobalVec); CHKERRQ(ierr);
ierr = DMPlexNaturalToGlobalEnd (ddm, lNatVec, lGlobalVec); CHKERRQ(ierr);
ierr = VecView(lNatVec, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
ierr = VecView(lGlobalVec, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
// Returning the Global vector to the natural one should give the initial values...
ierr = VecZeroEntries(lNatVec); CHKERRQ(ierr);
ierr = DMPlexGlobalToNaturalBegin(ddm, lGlobalVec, lNatVec); CHKERRQ(ierr);
ierr = DMPlexGlobalToNaturalEnd (ddm, lGlobalVec, lNatVec); CHKERRQ(ierr);
ierr = VecView(lNatVec, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
//Verification
ierr = VecGetArray(lNatVec, &lNatVecArray); CHKERRQ(ierr);
//Verify that we got the initial partition back into the "natural" vector:
PetscBool lResultBad = PETSC_FALSE;
for (c = 0; c < Nc; ++c) {
lResultBad |= lNatVecArray[c] != (InitPartForRank[rank])[c];
}
ierr = VecRestoreArray(lNatVec, &lNatVecArray); CHKERRQ(ierr);
// Now create an Overlap:
//ierr = DMSetBasicAdjacency(ddm, PETSC_TRUE, PETSC_FALSE); CHKERRQ(ierr); //This works: it give only 2 cells for the overlap
//ierr = DMSetAdjacency(ddm, 0, PETSC_TRUE, PETSC_FALSE); CHKERRQ(ierr); //This works too, but why do I have to use field #0 ??? Is it a convention?
ierr = DMSetAdjacency(ddm, PETSC_DEFAULT, PETSC_TRUE, PETSC_FALSE); CHKERRQ(ierr); //This does not works: it give 3 cells for the overlap instead of 2.
DM ddmo;
PetscSF sfo;
ierr = DMPlexDistributeOverlap(ddm, 1, &sfo, &ddmo); CHKERRQ(ierr);
IS lCellsWithOvl = 0;
ierr = DMPlexGetCellNumbering(ddmo, &lCellsWithOvl); CHKERRQ(ierr);
PetscScalar *lGlobVecArray = 0;
ierr = VecGetArray(lGlobalVec, &lGlobVecArray); CHKERRQ(ierr);
PetscInt n = 0;
const PetscInt *lCellsPetscNum = 0;
ierr = ISGetLocalSize(lCellsWithOvl, &n); CHKERRQ(ierr);
ierr = ISGetIndices(lCellsWithOvl, &lCellsPetscNum); CHKERRQ(ierr);
int i;
for (i = 0; i < n; i++) {
const int lGlobCellNum = lCellsPetscNum[i];
const int lLocalCellNum = lGlobCellNum - rank*5; //hard coded 5 elements by process
if (lGlobCellNum >= 0) {
const int lInitialCellNum = lGlobVecArray[lLocalCellNum];
printf("[%d] petsc_cell_# %d is initial (as in sInitialPartition) cell num %d\n", rank, lCellsPetscNum[i], lInitialCellNum);
}
else {
printf("[%d] overlap petsc_cell_# %d \n", rank, -lCellsPetscNum[i]-1);
}
}
ierr = ISRestoreIndices(lCellsWithOvl, &lCellsPetscNum);
ierr = ISView(lCellsWithOvl,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
ierr = VecRestoreArray(lGlobalVec, &lGlobVecArray); CHKERRQ(ierr);
ierr = DMRestoreGlobalVector(dm,&lNatVec); CHKERRQ(ierr);
ierr = DMRestoreGlobalVector(ddm,&lGlobalVec); CHKERRQ(ierr);
ierr = PetscSFDestroy(&sfMig);CHKERRQ(ierr);
ierr = PetscSFDestroy(&sfPart);CHKERRQ(ierr);
ierr = DMDestroy(&dm);CHKERRQ(ierr);
ierr = DMDestroy(&ddm);CHKERRQ(ierr);
ierr = PetscFinalize();
return ierr || lResultBad;
}
/*TEST
testset:
args: -field
nsize: 2
test:
suffix: 0
args:
test:
suffix: 1
args: -box -dm_plex_simplex 0 -dm_plex_box_faces 7,10 -dm_distribute
TEST*/
