Dear Koki,
Thanks for your email. In the example of your last email
DMPlexCoordinatesLoad() takes sF0 (PetscSF) as a third argument. In our
code this modification does not fix the error when loading a periodic
dm. Are we doing something wrong? I've included an example code at the
bottom of this email, including the error output.
Thanks and best regards,
Berend
/**** Write DM + Vec restart ****/
PetscViewerHDF5Open(PETSC_COMM_WORLD, "result", FILE_MODE_WRITE, &H5Viewer);
PetscObjectSetName((PetscObject)dm, "plexA");
PetscViewerPushFormat(H5Viewer, PETSC_VIEWER_HDF5_PETSC);
DMPlexTopologyView(dm, H5Viewer);
DMPlexLabelsView(dm, H5Viewer);
DMPlexCoordinatesView(dm, H5Viewer);
PetscViewerPopFormat(H5Viewer);
DM sdm;
PetscSection s;
DMClone(dm, &sdm);
PetscObjectSetName((PetscObject)sdm, "dmA");
DMGetGlobalSection(dm, &s);
DMSetGlobalSection(sdm, s);
DMPlexSectionView(dm, H5Viewer, sdm);
Vec vec, vecOld;
PetscScalar *array, *arrayOld, *xVecArray, *xVecArrayOld;
PetscInt numPoints;
DMGetGlobalVector(sdm, &vec);
DMGetGlobalVector(sdm, &vecOld);
/*** Fill the vectors vec and vecOld ***/
VecGetArray(vec, &array);
VecGetArray(vecOld, &arrayOld);
VecGetLocalSize(xGlobalVector, &numPoints);
VecGetArray(xGlobalVector, &xVecArray);
VecGetArray(xOldGlobalVector, &xVecArrayOld);
for (i = 0; i < numPoints; i++) /* Loop over all internal mesh points */
{
array[i] = xVecArray[i];
arrayOld[i] = xVecArrayOld[i];
}
VecRestoreArray(vec, &array);
VecRestoreArray(vecOld, &arrayOld);
VecRestoreArray(xGlobalVector, &xVecArray);
VecRestoreArray(xOldGlobalVector, &xVecArrayOld);
PetscObjectSetName((PetscObject)vec, "vecA");
PetscObjectSetName((PetscObject)vecOld, "vecB");
DMPlexGlobalVectorView(dm, H5Viewer, sdm, vec);
DMPlexGlobalVectorView(dm, H5Viewer, sdm, vecOld);
PetscViewerDestroy(&H5Viewer);
/*** end of writing ****/
/*** Load ***/
PetscViewerHDF5Open(PETSC_COMM_WORLD, "result", FILE_MODE_READ, &H5Viewer);
DMCreate(PETSC_COMM_WORLD, &dm);
DMSetType(dm, DMPLEX);
PetscObjectSetName((PetscObject)dm, "plexA");
PetscViewerPushFormat(H5Viewer, PETSC_VIEWER_HDF5_PETSC);
DMPlexTopologyLoad(dm, H5Viewer, &sfO);
DMPlexLabelsLoad(dm, H5Viewer);
DMPlexCoordinatesLoad(dm, H5Viewer, sfO);
PetscViewerPopFormat(H5Viewer);
DMPlexDistribute(dm, Options->Mesh.overlap, &sfDist, &distributedDM);
if (distributedDM) {
DMDestroy(&dm);
dm = distributedDM;
PetscObjectSetName((PetscObject)dm, "plexA");
}
PetscSFCompose(sfO, sfDist, &sf);
PetscSFDestroy(&sfO);
PetscSFDestroy(&sfDist);
DMClone(dm, &sdm);
PetscObjectSetName((PetscObject)sdm, "dmA");
DMPlexSectionLoad(dm, H5Viewer, sdm, sf, &globalDataSF, &localDataSF);
/** Load the Vectors **/
DMGetGlobalVector(sdm, &Restart_xGlobalVector);
VecSet(Restart_xGlobalVector,0.0);
PetscObjectSetName((PetscObject)Restart_xGlobalVector, "vecA");
DMPlexGlobalVectorLoad(dm, H5Viewer, sdm,
globalDataSF,Restart_xGlobalVector);
DMGetGlobalVector(sdm, &Restart_xOldGlobalVector);
VecSet(Restart_xOldGlobalVector,0.0);
PetscObjectSetName((PetscObject)Restart_xOldGlobalVector, "vecB");
DMPlexGlobalVectorLoad(dm, H5Viewer, sdm, globalDataSF,
Restart_xOldGlobalVector);
PetscViewerDestroy(&H5Viewer);
/**** The error message when loading is the following ************/
Creating and distributing mesh
[0]PETSC ERROR: --------------------- Error Message
--------------------------
[0]PETSC ERROR: Invalid argument
[0]PETSC ERROR: Number of coordinates loaded 17128 does not match number
of vertices 8000
[0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.
[0]PETSC ERROR: Petsc Development GIT revision: v3.16.1-435-g007f11b901
GIT Date: 2021-12-01 14:31:21 +0000
[0]PETSC ERROR: ./MF3 on a linux-gcc-openmpi-opt named
ivt24.ads.uni-magdeburg.de by berend Mon Dec 6 16:11:21 2021
[0]PETSC ERROR: Configure options --with-p4est=yes --with-partemis
--with-metis --with-debugging=no --download-metis=yes
--download-parmetis=yes --with-errorchecking=no --download-hdf5
--download-zlib --download-p4est
[0]PETSC ERROR: #1 DMPlexCoordinatesLoad_HDF5_V0_Private() at
/home/berend/src/petsc_main/src/dm/impls/plex/plexhdf5.c:1387
[0]PETSC ERROR: #2 DMPlexCoordinatesLoad_HDF5_Internal() at
/home/berend/src/petsc_main/src/dm/impls/plex/plexhdf5.c:1419
[0]PETSC ERROR: #3 DMPlexCoordinatesLoad() at
/home/berend/src/petsc_main/src/dm/impls/plex/plex.c:2070
[0]PETSC ERROR: #4 RestartMeshDM() at
/home/berend/src/eclipseworkspace/multiflow/src/io/restartmesh.c:81
[0]PETSC ERROR: #5 CreateMeshDM() at
/home/berend/src/eclipseworkspace/multiflow/src/mesh/createmesh.c:61
[0]PETSC ERROR: #6 main() at
/home/berend/src/eclipseworkspace/multiflow/src/general/main.c:132
[0]PETSC ERROR: PETSc Option Table entries:
[0]PETSC ERROR: --download-hdf5
[0]PETSC ERROR: --download-metis=yes
[0]PETSC ERROR: --download-p4est
[0]PETSC ERROR: --download-parmetis=yes
[0]PETSC ERROR: --download-zlib
[0]PETSC ERROR: --with-debugging=no
[0]PETSC ERROR: --with-errorchecking=no
[0]PETSC ERROR: --with-metis
[0]PETSC ERROR: --with-p4est=yes
[0]PETSC ERROR: --with-partemis
[0]PETSC ERROR: -d results
[0]PETSC ERROR: -o run.mf
[0]PETSC ERROR: ----------------End of Error Message -------send entire
error message to [email protected]
--------------------------------------------------------------------------
MPI_ABORT was invoked on rank 0 in communicator MPI_COMM_WORLD
with errorcode 62.
NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
You may or may not see output from other processes, depending on
exactly when Open MPI kills them.
--------------------------------------------------------------------------
On 11/19/21 00:26, Sagiyama, Koki wrote:
Hi Berend,
I was not able to reproduce the issue you are having, but the following
1D example (and similar 2D examples) worked fine for me using the latest
PETSc. Please note that DMPlexCoordinatesLoad() now takes a PetscSF
object as the third argument, but the default behavior is unchanged.
/* test_periodic_io.c */
#include <petscdmplex.h>
#include <petscsf.h>
#include <petsclayouthdf5.h>
int main(int argc, char **argv)
{
DM dm;
Vec coordinates;
PetscViewer viewer;
PetscViewerFormat format = PETSC_VIEWER_HDF5_PETSC;
PetscSF sfO;
PetscErrorCode ierr;
ierr = PetscInitialize(&argc, &argv, NULL, NULL); if (ierr) return ierr;
/* Save */
ierr = PetscViewerHDF5Open(PETSC_COMM_WORLD, "periodic_example.h5",
FILE_MODE_WRITE, &viewer);CHKERRQ(ierr);
{
DM pdm;
PetscInt dim = 1;
const PetscInt faces[1] = {4};
DMBoundaryType periodicity[] = {DM_BOUNDARY_PERIODIC};
PetscInt overlap = 1;
ierr = DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, PETSC_FALSE,
faces, NULL, NULL, periodicity, PETSC_TRUE, &dm);CHKERRQ(ierr);
ierr = DMPlexDistribute(dm, overlap, NULL, &pdm);CHKERRQ(ierr);
ierr = DMDestroy(&dm);CHKERRQ(ierr);
dm = pdm;
ierr = PetscObjectSetName((PetscObject)dm, "periodicDM");CHKERRQ(ierr);
}
ierr = DMGetCoordinates(dm, &coordinates);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD, "Coordinates before
saving:\n");CHKERRQ(ierr);
ierr = VecView(coordinates, NULL);CHKERRQ(ierr);
ierr = PetscViewerPushFormat(viewer, format);CHKERRQ(ierr);
ierr = DMPlexTopologyView(dm, viewer);CHKERRQ(ierr);
ierr = DMPlexCoordinatesView(dm, viewer);CHKERRQ(ierr);
ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
ierr = DMDestroy(&dm);CHKERRQ(ierr);
ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
/* Load */
ierr = PetscViewerHDF5Open(PETSC_COMM_WORLD, "periodic_example.h5",
FILE_MODE_READ, &viewer);CHKERRQ(ierr);
ierr = DMCreate(PETSC_COMM_WORLD, &dm);CHKERRQ(ierr);
ierr = DMSetType(dm, DMPLEX);CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject)dm, "periodicDM");CHKERRQ(ierr);
ierr = PetscViewerPushFormat(viewer, format);CHKERRQ(ierr);
ierr = DMPlexTopologyLoad(dm, viewer, &sfO);CHKERRQ(ierr);
ierr = DMPlexCoordinatesLoad(dm, viewer, sfO);CHKERRQ(ierr);
ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
ierr = DMGetCoordinates(dm, &coordinates);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD, "Coordinates after
loading:\n");CHKERRQ(ierr);
ierr = VecView(coordinates, NULL);CHKERRQ(ierr);
ierr = PetscSFDestroy(&sfO);CHKERRQ(ierr);
ierr = DMDestroy(&dm);CHKERRQ(ierr);
ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
ierr = PetscFinalize();
return ierr;
}
mpiexec -n 2 ./test_periodic_io
Coordinates before saving:
Vec Object: coordinates 2 MPI processes
type: mpi
Process [0]
0.
Process [1]
0.25
0.5
0.75
Coordinates after loading:
Vec Object: vertices 2 MPI processes
type: mpi
Process [0]
0.
0.25
0.5
0.75
Process [1]
I would also like to note that, with the latest update, we can
optionally load coordinates directly on the distributed dm as (using
your notation):
/* Distribute dm */
...
PetscSFCompose(sfO, sfDist, &sf);
DMPlexCoordinatesLoad(dm, viewer, sf);
To use this feature, we need to pass "-dm_plex_view_hdf5_storage_version
2.0.0" option when saving topology/coordinates.
Thanks,
Koki
------------------------------------------------------------------------
*From:* Berend van Wachem <[email protected]>
*Sent:* Wednesday, November 17, 2021 3:16 PM
*To:* Hapla Vaclav <[email protected]>; PETSc users list
<[email protected]>; Lawrence Mitchell <[email protected]>; Sagiyama,
Koki <[email protected]>
*Subject:* Re: [petsc-users] DMView and DMLoad
*******************
This email originates from outside Imperial. Do not click on links and
attachments unless you recognise the sender.
If you trust the sender, add them to your safe senders list
https://spam.ic.ac.uk/SpamConsole/Senders.aspx
<https://spam.ic.ac.uk/SpamConsole/Senders.aspx> to disable email
stamping for this address.
*******************
Dear Vaclav, Lawrence, Koki,
Thanks for your help! Following your advice and following your example
(https://petsc.org/main/docs/manual/dmplex/#saving-and-loading-data-with-hdf5
<https://petsc.org/main/docs/manual/dmplex/#saving-and-loading-data-with-hdf5>)
we are able to save and load the DM with a wrapped Vector in h5 format
(PETSC_VIEWER_HDF5_PETSC) successfully.
For saving, we use something similar to:
DMPlexTopologyView(dm, viewer);
DMClone(dm, &sdm);
...
DMPlexSectionView(dm, viewer, sdm);
DMGetLocalVector(sdm, &vec);
...
DMPlexLocalVectorView(dm, viewer, sdm, vec);
and for loading:
DMCreate(PETSC_COMM_WORLD, &dm);
DMSetType(dm, DMPLEX);
...
PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_PETSC);
DMPlexTopologyLoad(dm, viewer, &sfO);
DMPlexLabelsLoad(dm, viewer);
DMPlexCoordinatesLoad(dm, viewer);
PetscViewerPopFormat(viewer);
...
PetscSFCompose(sfO, sfDist, &sf);
...
DMClone(dm, &sdm);
DMPlexSectionLoad(dm, viewer, sdm, sf, &globalDataSF, &localDataSF);
DMGetLocalVector(sdm, &vec);
...
DMPlexLocalVectorLoad(dm, viewer, sdm, localDataSF, vec);
This works fine for non-periodic DMs but for periodic cases the line:
DMPlexCoordinatesLoad(dm, H5Viewer);
delivers the error message: invalid argument and the number of loaded
coordinates does not match the number of vertices.
Is this a known shortcoming, or have we forgotten something to load
periodic DMs?
Best regards,
Berend.
On 9/22/21 20:59, Hapla Vaclav wrote:
To avoid confusions here, Berend seems to be specifically demanding XDMF
(PETSC_VIEWER_HDF5_XDMF). The stuff we are now working on is parallel
checkpointing in our own HDF5 format (PETSC_VIEWER_HDF5_PETSC), I will
make a series of MRs on this topic in the following days.
For XDMF, we are specifically missing the ability to write/load DMLabels
properly. XDMF uses specific cell-local numbering for faces for
specification of face sets, and face-local numbering for specification
of edge sets, which is not great wrt DMPlex design. And ParaView doesn't
show any of these properly so it's hard to debug. Matt, we should talk
about this soon.
Berend, for now, could you just load the mesh initially from XDMF and
then use our PETSC_VIEWER_HDF5_PETSC format for subsequent saving/loading?
Thanks,
Vaclav
On 17 Sep 2021, at 15:46, Lawrence Mitchell <[email protected]
<mailto:[email protected] <mailto:[email protected]>>> wrote:
Hi Berend,
On 14 Sep 2021, at 12:23, Matthew Knepley <[email protected]
<mailto:[email protected] <mailto:[email protected]>>> wrote:
On Tue, Sep 14, 2021 at 5:15 AM Berend van Wachem
<[email protected] <mailto:[email protected] <mailto:[email protected]>>> wrote:
Dear PETSc-team,
We are trying to save and load distributed DMPlex and its associated
physical fields (created with DMCreateGlobalVector) (Uvelocity,
VVelocity, ...) in HDF5_XDMF format. To achieve this, we do the
following:
1) save in the same xdmf.h5 file:
DMView( DM , H5_XDMF_Viewer );
VecView( UVelocity, H5_XDMF_Viewer );
2) load the dm:
DMPlexCreateFromfile(PETSC_COMM_WORLD, Filename, PETSC_TRUE, DM);
3) load the physical field:
VecLoad( UVelocity, H5_XDMF_Viewer );
There are no errors in the execution, but the loaded DM is distributed
differently to the original one, which results in the incorrect
placement of the values of the physical fields (UVelocity etc.) in the
domain.
This approach is used to restart the simulation with the last saved DM.
Is there something we are missing, or there exists alternative routes to
this goal? Can we somehow get the IS of the redistribution, so we can
re-distribute the vector data as well?
Many thanks, best regards,
Hi Berend,
We are in the midst of rewriting this. We want to support saving
multiple meshes, with fields attached to each,
and preserving the discretization (section) information, and allowing
us to load up on a different number of
processes. We plan to be done by October. Vaclav and I are doing this
in collaboration with Koki Sagiyama,
David Ham, and Lawrence Mitchell from the Firedrake team.
The core load/save cycle functionality is now in PETSc main. So if
you're using main rather than a release, you can get access to it now.
This section of the manual shows an example of how to do
thingshttps://petsc.org/main/docs/manual/dmplex/#saving-and-loading-data-with-hdf5
<https://petsc.org/main/docs/manual/dmplex/#saving-and-loading-data-with-hdf5
<https://petsc.org/main/docs/manual/dmplex/#saving-and-loading-data-with-hdf5>>
Let us know if things aren't clear!
Thanks,
Lawrence
