[petsc-users] Write/read a DMForest to/from a file

[Berend van Wachem]

Dear Mark,

Yes, I have tried that. That will only work if I convert the DMForest to 
a DMPlex first, and then write the DMPlex to file. But then, I cannot 
read in a DMForest when I want to continue the calculations later on.

Below is the code I use to write a DMPlex to file. When I call this 
routine with a DMForest, I get the error:

[0]PETSC ERROR: Invalid argument
[0]PETSC ERROR: Wrong type of object: Parameter # 1
[0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.18.2, Nov 28, 2022
[0]PETSC ERROR: ./bin/write_periodic on a linux-gcc-dev named 
multiflow.multiflow.org by berend Fri Dec 30 11:54:51 2022
[0]PETSC ERROR: Configure options --with-clean --download-metis=yes 
--download-parmetis=yes --download-hdf5 --download-p4est 
--download-triangle --download-tetgen --with-zlib-lib=/usr/lib64/libz.a 
--with-zlib-include=/usr/include --with-mpi=yes --with-mpi-dir=/usr 
--with-mpiexec=/usr/bin/mpiexec --download-slepc=yes --download-fftw=yes
[0]PETSC ERROR: #1 PetscSectionGetChart() at 
/usr/local/petsc-3.18.2/src/vec/is/section/interface/section.c:592
[0]PETSC ERROR: #2 DMPlexGetChart() at 
/usr/local/petsc-3.18.2/src/dm/impls/plex/plex.c:2702
[0]PETSC ERROR: #3 DMPlexGetDepthStratum() at 
/usr/local/petsc-3.18.2/src/dm/impls/plex/plex.c:4882
[0]PETSC ERROR: #4 DMPlexCreatePointNumbering() at 
/usr/local/petsc-3.18.2/src/dm/impls/plex/plex.c:8303
[0]PETSC ERROR: #5 DMPlexTopologyView() at 
/usr/local/petsc-3.18.2/src/dm/impls/plex/plex.c:1854
[0]PETSC ERROR: #6 WriteDMAndVectortoHDF5File() at 
/home/berend/git/Code/petsc-dmplex-restart-test/src/createandwrite.c:176
[0]PETSC ERROR: #7 main() at 
/home/berend/git/Code/petsc-dmplex-restart-test/src/createandwrite.c:555

The code I use to write is:

PetscErrorCode WriteDMAndVectortoHDF5File(DM dm, Vec *DataVector, const 
char *HDF5file)
{
   PetscViewer H5Viewer;
   PetscErrorCode ierr;
   DM dmCopy;
   PetscSection sectionCopy;
   PetscInt i;
   PetscScalar *xVecArray;
   PetscInt numPoints, numPointsCopy;
   Vec vectorCopy;
   PetscScalar *array;

   PetscFunctionBegin;
   ierr = PetscViewerHDF5Open(PETSC_COMM_WORLD, HDF5file, 
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);

   ierr = DMClone(dm, &dmCopy);
   CHKERRQ(ierr);
   ierr = PetscObjectSetName((PetscObject) dmCopy, "DM");
   CHKERRQ(ierr);
   ierr = DMGetLocalSection(dm, &sectionCopy);
   CHKERRQ(ierr);
   ierr = DMSetLocalSection(dmCopy, sectionCopy);
   CHKERRQ(ierr);

   /* Write the section to the file */
   ierr = DMPlexSectionView(dm, H5Viewer, dmCopy);
   CHKERRQ(ierr);

   ierr = DMGetGlobalVector(dmCopy, &vectorCopy);
   CHKERRQ(ierr);

   /*** We have to copy the vector into the new vector ... ***/
   ierr = VecGetArray(vectorCopy, &array);
   CHKERRQ(ierr);
   ierr = VecGetLocalSize(*DataVector, &numPoints);
   CHKERRQ(ierr);
   ierr = VecGetLocalSize(vectorCopy, &numPointsCopy);
   CHKERRQ(ierr);
   assert(numPoints == numPointsCopy);
   ierr = VecGetArray(*DataVector, &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(*DataVector, &xVecArray);
   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);

   ierr = DMDestroy(&dmCopy);
   CHKERRQ(ierr);
   /*** End of writing ****/
   PetscFunctionReturn(0);
}





On 29/12/2022 21:54, Mark Adams wrote:
> Have you tried using the DMForest as you would use a DMPLex?
> That should work. Forst has extra stuff but it just makes a Plex and the 
> end of the day.
> 
> On Tue, Dec 27, 2022 at 5:21 AM Berend van Wachem 
> <berend.vanwachem at ovgu.de <mailto:berend.vanwachem at ovgu.de>> wrote:
> 
>     Dear Petsc-team/users,
> 
>     I am trying to save a DMForest which has been used for a calculation
>     (so
>     refined/coarsened in places) to a file and later read it from the file
>     to continue a calculation with.
>     My question is: how do I do this? I've tried a few things, such as
>     using
>     DMView and DMLoad, and saving the BaseDM, but that doesn't seem to save
>     the refining/coarsening which has taken place in the initial
>     calculation.
>     I haven't been able to find any example or documentation for this. Any
>     pointers or examples would be very much appreciated!
> 
>     Best regards, Berend.
>