[petsc-users] VecLoad from HDF5 file
Matthew Knepley
knepley at gmail.com
Thu Apr 22 16:40:43 CDT 2021
On Thu, Apr 22, 2021 at 3:00 PM Hapla Vaclav <vaclav.hapla at erdw.ethz.ch>
wrote:
> Thibault, would it make sense to you if we replace the current
> DMSetOutputSequenceNumber() and PetscViewerHDF5SetTimestep() by just
> VecSetTimestep()? Details here
> <https://gitlab.com/petsc/petsc/-/merge_requests/3881#note_558257800>.
>
> From what I can see, this could simplify and robustify the code and reduce
> docs scattering.
>
> Matt, it seems that the OutputSequenceNumber was meant as some more
> general concept, but now it's exclusively interpreted as a timestep
> selector for Vec I/O. See
> grep -RIn -e DMGetOutputSequenceNumber -e DMSetOutputSequenceNumber -e
> outputSequenceNum src | less
>
Well, I was thinking that we would need sequencing for things not currently
on our agenda, such as optimization steps, or design iterations, etc. so I
wanted to keep the terminology somewhat open.
Thanks,
Matt
> Thanks,
> Vaclav
>
> On 9 Apr 2021, at 16:34, Thibault Bridel-Bertomeu <
> thibault.bridelbertomeu at gmail.com> wrote:
>
> Yes hmmm more documentation could simply be added to VecLoad ? What do you
> think ? I mean there are already a few special cases mentioned in there, so
> maybe adding the case where the Vec stems from a DM could be helpful.
> What's more, I think VecLoad is the first function users are going to look
> up before going into more details, so it could as well have as much details
> as possible.
>
> Anyways if anyone looks up this thread in the mailing list, here are the
> snippets that do the job.
>
> First snippet : how to write a HDF5 file with only the last iteration of a
> solution vector carried by a DM and a TS :
>
> call DMGetOutputSequenceNumber(dm, save_seqnum, save_seqval, ierr);
> CHKERRA(ierr)
> call DMSetOutputSequenceNumber(dm, -1, 0.d0, ierr); CHKERRA(ierr)
> call PetscViewerCreate(PETSC_COMM_WORLD, hdf5Viewer, ierr); CHKERRA(ierr)
> call PetscViewerSetType(hdf5Viewer, PETSCVIEWERHDF5, ierr); CHKERRA(ierr)
> call PetscViewerFileSetMode(hdf5Viewer, FILE_MODE_WRITE, ierr);
> CHKERRA(ierr);
> write(filename,'(A,I5.5,A)') "restart_", stepnum, ".h5"
> call PetscViewerFileSetName(hdf5Viewer, trim(filename), ierr);
> CHKERRA(ierr)
> call VecView(X, hdf5Viewer, ierr); CHKERRA(ierr)
> call PetscViewerDestroy(hdf5Viewer, ierr); CHKERRA(ierr)
> call DMSetOutputSequenceNumber(dm, save_seqnum, save_seqval, ierr);
> CHKERRA(ierr)
>
> Second snippet : how to read such a HDF5 file back into a new DM that has
> no idea whatsoever what happened when the file was written :
>
> call DMGetOutputSequenceNumber(dm, save_seqnum, save_seqval, ierr);
> CHKERRA(ierr)
> call DMSetOutputSequenceNumber(dm, -1, 0.d0, ierr); CHKERRA(ierr)
> call PetscViewerHDF5Open(PETSC_COMM_WORLD, trim(restartname),
> FILE_MODE_READ, hdf5Viewer, ierr); CHKERRA(ierr)
> call PetscViewerHDF5SetTimestep(hdf5Viewer, restartiter, ierr);
> CHKERRA(ierr)
> call VecLoad(sol, hdf5Viewer, ierr); CHKERRA(ierr)
> call PetscViewerDestroy(hdf5Viewer, ierr); CHKERRA(ierr)
> call DMSetOutputSequenceNumber(dm, save_seqnum, save_seqval, ierr);
> CHKERRA(ierr)
>
> Note that in both snippets, vectors X and sol resp. are global vectors.
>
> Cheers and thanks again for the help !!
>
> Thibault
>
> Le jeu. 8 avr. 2021 à 22:17, Matthew Knepley <knepley at gmail.com> a écrit :
>
>> On Thu, Apr 8, 2021 at 9:00 AM Thibault Bridel-Bertomeu <
>> thibault.bridelbertomeu at gmail.com> wrote:
>>
>>> Hello everyone,
>>>
>>> So i figured it out, and I got a work-around, but it's not pretty and I
>>> think there might be a nicer solution that you guys will see.
>>>
>>> Here is the stack of calls for the VecLoad of a vector generated via a
>>> DMPlex :
>>>
>>> CALLING PETSCVIEWERHDF5SETTIMESTEP FROM FORTRAN
>>> TOTOTOTOTO :::>>> IN PETSCVIEWERHDF5SETTIMESTEP. <hdf5v.c>
>>> TOTOTOTOTO :::>>> HDF5->TIMESTEP = 20. <hdf5v.c>
>>> CALLING VECLOAD FROM FORTRAN
>>> TOTOTOTOTO :::>>> INSIDE VECLOAD <vector.c>
>>> TOTOTOTOTO :::>>> FORMAT = PETSC_VIEWER_RAW. <vector.c>
>>> TOTOTOTOTO :::>>> IN VECLOAD_PLEX. <plex.c>
>>> TOTOTOTOTO :::>>> IN VECLOAD_PLEX_HDF5_INTERNAL <plexhdf5.c>
>>> TOTOTOTOTO :::>>> IN PETSCVIEWERHDF5SETTIMESTEP. <hdf5v.c>
>>> TOTOTOTOTO :::>>> HDF5->TIMESTEP = 20. <hdf5v.c>
>>> TOTOTOTOTO :::>>> IN VECLOAD_DEFAULT. <vecio.c>
>>> TOTOTOTOTO :::>>> IN VECLOAD_HDF5 <vecio.c>.
>>> TOTOTOTOTO :::>>> IN PETSCVIEWERHDF5LOAD. <hdf5io.c>
>>> TOTOTOTOTO :::>>> IN PETSCVIEWERHDF5READINITIALIZE_PRIVATE. <hdf5io.c>
>>> TOTOTOTOTO :::>>> TIMESTEP = 20. <hdf5io.c>
>>> TOTOTOTOTO :::>>> LENIND = 1, TIMESTEP = 20. <hdf5io.c>
>>>
>>> Not everything is there, I just stopped at line 78 of hdf5io.c as Matt
>>> suggested.
>>> The thing is, you can call
>>>
>>> call PetscViewerHDF5SetTimestep(hdf5Viewer, restartiter, ierr);
>>> CHKERRA(ierr)
>>>
>>> all you want before calling
>>>
>>> call VecLoad(sol, hdf5Viewer, ierr); CHKERRA(ierr)
>>>
>>> it won't have any effect.
>>> The timestep that is set by the ..SetTimestep routine is overriden
>>> inside the
>>>
>>> VecLoad_Plex_HDF5_Internal
>>>
>>> routine because of these lines :
>>>
>>> ierr = DMGetOutputSequenceNumber(dm, &seqnum, NULL);CHKERRQ(ierr);
>>> ierr = PetscViewerHDF5SetTimestep(viewer, seqnum);CHKERRQ(ierr);
>>>
>>> Since my point is to "restart" a computation from such a HDF5 file, when
>>> i read the file the DMPlex is brand new and has no idea what happened
>>> during another run of the program --> the "seqnum" is set to -1.
>>> As a consequence, the viewer always sees its "timestep" variable as
>>> being -1 even if you try forcing it with
>>>
>>> call PetscViewerHDF5SetTimestep(hdf5Viewer, restartiter, ierr);
>>> CHKERRA(ierr)
>>>
>>> Anyways, one way to do it is then to add :
>>>
>>> call DMSetOutputSequenceNumber(dm, restartiter, PETSC_NULL_REAL, ierr);
>>> CHKERRA(ierr)
>>>
>>> to the code, that thus becomes :
>>>
>>> call DMSetOutputSequenceNumber(dm, restartiter, PETSC_NULL_REAL, ierr);
>>> CHKERRA(ierr)
>>> call PetscViewerHDF5Open(PETSC_COMM_WORLD, trim(restartname),
>>> FILE_MODE_READ, hdf5Viewer, ierr); CHKERRA(ierr)
>>> call PetscViewerHDF5PushGroup(hdf5Viewer, "/fields", ierr);
>>> CHKERRA(ierr)
>>> call PetscViewerHDF5SetTimestep(hdf5Viewer, restartiter, ierr);
>>> CHKERRA(ierr)
>>> call VecLoad(sol, hdf5Viewer, ierr); CHKERRA(ierr)
>>> call PetscViewerDestroy(hdf5Viewer, ierr); CHKERRA(ierr)
>>>
>>> By the way, it even works if you remove the PushGroup(..., "/fields",
>>> ...) because this is also forced in VecLoad_Plex_HDF5_Internal.
>>>
>>> Might be a better idea would be to do a DMLoad, what do you think ?
>>>
>>
>> Okay, this is correct. I should have thought of this. If a Vec comes from
>> a DM, then the DM controls the I/O. We should be
>> more explicit about this in the documentation. From this decision flows
>> the decision to have the Solver (TS) communicate
>> with the DM instead of the Viewer.
>>
>> Thus, if you want to load a certain output in the absence of a solver,
>> you too should communicate with it through the DM. I think
>> the problem here is lack of documentation rather than the choice. I will
>> think about where to document this. If you have a suggestion,
>> maybe where you looked first, it would be appreciated.
>>
>> Thanks,
>>
>> Matt
>>
>>
>>> Thanks !!
>>>
>>> Thibault
>>>
>>> Le jeu. 8 avr. 2021 à 13:35, Thibault Bridel-Bertomeu <
>>> thibault.bridelbertomeu at gmail.com> a écrit :
>>>
>>>> Argh no I don't, I implemented it straight in a larger CFD code because
>>>> everything is already there for setting-up the DM, DS, PetscFV etc...
>>>> Now that you guys gave me some more pointers I'll dig around with a
>>>> debugger ..
>>>>
>>>> Thanks !
>>>>
>>>> Le jeu. 8 avr. 2021 à 13:25, Matthew Knepley <knepley at gmail.com> a
>>>> écrit :
>>>>
>>>>> Do you have a self-contained thing I could run? This should be easy to
>>>>> clear up in the debugger.
>>>>>
>>>>> Thanks,
>>>>>
>>>>> Matt
>>>>>
>>>>> On Thu, Apr 8, 2021 at 7:22 AM Thibault Bridel-Bertomeu <
>>>>> thibault.bridelbertomeu at gmail.com> wrote:
>>>>>
>>>>>> Doesn't seem to be a Fortran stub problem, it's automatically
>>>>>> generated and looks good to me ...
>>>>>>
>>>>>>
>>>>>>
>>>>>> Le jeu. 8 avr. 2021 à 13:20, Thibault Bridel-Bertomeu <
>>>>>> thibault.bridelbertomeu at gmail.com> a écrit :
>>>>>>
>>>>>>> Well ... I cleaned everything and did a fresh compile just now with
>>>>>>> "print*, "foo"" here and there and it does print foo on the terminal ...
>>>>>>> So yes, I am running the changed code ! ;)
>>>>>>>
>>>>>>> call DMCreateGlobalVector(dm, sol, ierr); CHKERRA(ierr)
>>>>>>> call VecZeroEntries(sol, ierr); CHKERRA(ierr)
>>>>>>> call PetscObjectSetName(sol, "Solution", ierr); CHKERRA(ierr)
>>>>>>> print*, "foo"
>>>>>>> call PetscViewerCreate(PETSC_COMM_WORLD, hdf5Viewer, ierr);
>>>>>>> CHKERRA(ierr)
>>>>>>> call PetscViewerSetType(hdf5Viewer, PETSCVIEWERHDF5, ierr);
>>>>>>> CHKERRA(ierr)
>>>>>>> call PetscViewerFileSetMode(hdf5Viewer, FILE_MODE_READ, ierr);
>>>>>>> CHKERRA(ierr)
>>>>>>> call PetscViewerFileSetName(hdf5Viewer, trim(restartname), ierr);
>>>>>>> CHKERRA(ierr)
>>>>>>> call PetscViewerHDF5SetTimestep(hdf5Viewer, restartiter, ierr);
>>>>>>> CHKERRA(ierr)
>>>>>>> call PetscViewerHDF5PushGroup(hdf5Viewer, "/fields", ierr);
>>>>>>> CHKERRA(ierr)
>>>>>>> call VecLoad(sol, hdf5Viewer, ierr); CHKERRA(ierr)
>>>>>>> call PetscViewerDestroy(hdf5Viewer, ierr); CHKERRA(ierr)
>>>>>>>
>>>>>>> In the snippet, restartiter is an integer equal to 20.
>>>>>>> I agree with you the line you just sent me should fix the issue.
>>>>>>> Maybe it's a fortran stub problem, I'll double check the stubs.
>>>>>>>
>>>>>>> Thanks !
>>>>>>> Thibault
>>>>>>>
>>>>>>> Le jeu. 8 avr. 2021 à 13:11, Matthew Knepley <knepley at gmail.com> a
>>>>>>> écrit :
>>>>>>>
>>>>>>>> On Thu, Apr 8, 2021 at 4:08 AM Thibault Bridel-Bertomeu <
>>>>>>>> thibault.bridelbertomeu at gmail.com> wrote:
>>>>>>>>
>>>>>>>>> Hi everyone,
>>>>>>>>>
>>>>>>>>> Thank you for your answers.
>>>>>>>>> Unfortunately, it does not work yet.
>>>>>>>>>
>>>>>>>>> 1/ I first tried just adding
>>>>>>>>>
>>>>>>>>> call PetscViewerHDF5SetTimestep(hdf5Viewer, 20, ierr);
>>>>>>>>> CHKERRA(ierr)
>>>>>>>>>
>>>>>>>>> (the latest iteration in the *.hdf5 is iter 20) before the
>>>>>>>>> VecLoad, i.e.:
>>>>>>>>>
>>>>>>>>> call PetscViewerCreate(PETSC_COMM_WORLD, hdf5Viewer, ierr);
>>>>>>>>> CHKERRA(ierr)
>>>>>>>>> call PetscViewerSetType(hdf5Viewer, PETSCVIEWERHDF5, ierr);
>>>>>>>>> CHKERRA(ierr)
>>>>>>>>> call PetscViewerFileSetMode(hdf5Viewer, FILE_MODE_READ, ierr);
>>>>>>>>> CHKERRA(ierr)
>>>>>>>>> call PetscViewerFileSetName(hdf5Viewer, trim(restartname), ierr);
>>>>>>>>> CHKERRA(ierr)
>>>>>>>>> call PetscViewerHDF5SetTimestep(hdf5Viewer, 20, ierr);
>>>>>>>>> CHKERRA(ierr)
>>>>>>>>> call PetscViewerHDF5PushGroup(hdf5Viewer, "/fields", ierr);
>>>>>>>>> CHKERRA(ierr)
>>>>>>>>> call VecLoad(sol, hdf5Viewer, ierr); CHKERRA(ierr)
>>>>>>>>> call PetscViewerDestroy(hdf5Viewer, ierr); CHKERRA(ierr)
>>>>>>>>>
>>>>>>>>> but it does not change anything, the error is exactly the same as
>>>>>>>>> before (with or without the PushGroup by the way).
>>>>>>>>>
>>>>>>>>
>>>>>>>> This is the one I do not believe. If the Viewer has a non-negative
>>>>>>>> timestep, it will increase the length index:
>>>>>>>>
>>>>>>>>
>>>>>>>> https://gitlab.com/petsc/petsc/-/blob/main/src/vec/is/utils/hdf5io.c#L78
>>>>>>>>
>>>>>>>> and we can see all the sizes in your HDF5 file. Are you sure you
>>>>>>>> are running the one with the changed code?
>>>>>>>>
>>>>>>>> Thanks,
>>>>>>>>
>>>>>>>> Matt
>>>>>>>>
>>>>>>>>
>>>>>>>>>
>>>>>>>>> 2/ I then tried adding
>>>>>>>>>
>>>>>>>>> call PetscViewerHDF5SetTimestep(hdf5Viewer, -1, ierr);
>>>>>>>>> CHKERRA(ierr)
>>>>>>>>>
>>>>>>>>> before the VecView to "prevent blocking with timesteps" as the doc
>>>>>>>>> says (
>>>>>>>>> https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Viewer/PetscViewerHDF5SetTimestep.html#PetscViewerHDF5SetTimestep),
>>>>>>>>> i.e. :
>>>>>>>>>
>>>>>>>>> call PetscViewerCreate(PETSC_COMM_WORLD, hdf5Viewer, ierr);
>>>>>>>>> CHKERRA(ierr)
>>>>>>>>> call PetscViewerSetType(hdf5Viewer, PETSCVIEWERHDF5, ierr);
>>>>>>>>> CHKERRA(ierr)
>>>>>>>>> call PetscViewerFileSetMode(hdf5Viewer, FILE_MODE_WRITE, ierr);
>>>>>>>>> CHKERRA(ierr);
>>>>>>>>> write(filename,'(A,I5.5,A)') "restart_", stepnum, ".h5"
>>>>>>>>> call PetscViewerFileSetName(hdf5Viewer, trim(filename), ierr);
>>>>>>>>> CHKERRA(ierr)
>>>>>>>>> call PetscViewerHDF5SetTimestep(hdf5Viewer, -1, ierr);
>>>>>>>>> CHKERRA(ierr)
>>>>>>>>> call VecView(X, hdf5Viewer, ierr); CHKERRA(ierr)
>>>>>>>>> call PetscViewerDestroy(hdf5Viewer, ierr); CHKERRA(ierr)
>>>>>>>>>
>>>>>>>>> but the hdf5 files keep increasing in size and each one stores
>>>>>>>>> more data than the previous one ... can I not disable that "stacking"
>>>>>>>>> effect and just keep the iteration I want ?
>>>>>>>>>
>>>>>>>>> 3/ I tried adding
>>>>>>>>>
>>>>>>>>> call PetscViewerHDF5SetTimestep(hdf5Viewer, stepnum, ierr);
>>>>>>>>> CHKERRA(ierr)
>>>>>>>>>
>>>>>>>>> before the VecView to try and just keep one iteration in the file,
>>>>>>>>> i.e. :
>>>>>>>>>
>>>>>>>>> call PetscViewerCreate(PETSC_COMM_WORLD, hdf5Viewer, ierr);
>>>>>>>>> CHKERRA(ierr)
>>>>>>>>> call PetscViewerSetType(hdf5Viewer, PETSCVIEWERHDF5, ierr);
>>>>>>>>> CHKERRA(ierr)
>>>>>>>>> call PetscViewerFileSetMode(hdf5Viewer, FILE_MODE_WRITE, ierr);
>>>>>>>>> CHKERRA(ierr);
>>>>>>>>> write(filename,'(A,I5.5,A)') "restart_", stepnum, ".h5"
>>>>>>>>> call PetscViewerFileSetName(hdf5Viewer, trim(filename), ierr);
>>>>>>>>> CHKERRA(ierr)
>>>>>>>>> call PetscViewerHDF5SetTimestep(hdf5Viewer, stepnum, ierr);
>>>>>>>>> CHKERRA(ierr)
>>>>>>>>> call VecView(X, hdf5Viewer, ierr); CHKERRA(ierr)
>>>>>>>>> call PetscViewerDestroy(hdf5Viewer, ierr); CHKERRA(ierr)
>>>>>>>>>
>>>>>>>>> but it has no effect whatsoever on the hdf5 files that are written
>>>>>>>>> ... they still store several timesteps, and then I cannot VecLoad what I
>>>>>>>>> want, i still get the same error.
>>>>>>>>>
>>>>>>>>> Any more idea ? I keep reading the doc on the HDF5 viewer and
>>>>>>>>> thinking that this or that function could solve the issue but either it has
>>>>>>>>> zero effect or does not solve the issue ...
>>>>>>>>>
>>>>>>>>> Thanks !!!
>>>>>>>>>
>>>>>>>>> Thibault
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> Le jeu. 8 avr. 2021 à 02:42, Barry Smith <bsmith at petsc.dev> a
>>>>>>>>> écrit :
>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> Matt,
>>>>>>>>>>
>>>>>>>>>> Is there anyway to provide this "extra meta-data" inside the
>>>>>>>>>> generated HDF file? Then when a fresh viewer opens the file it uses this
>>>>>>>>>> meta-data to know that the file contains "time-steps" and allows processing
>>>>>>>>>> of them? The simplest thing would be to have the viewer generate an error
>>>>>>>>>> if the file has time-steps but the user does not "request a (or more)
>>>>>>>>>> time-steps" instead of getting confused about the vector sizes.
>>>>>>>>>>
>>>>>>>>>> Barry
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> On Apr 7, 2021, at 6:32 PM, Matthew Knepley <knepley at gmail.com>
>>>>>>>>>> wrote:
>>>>>>>>>>
>>>>>>>>>> On Wed, Apr 7, 2021 at 4:10 PM Thibault Bridel-Bertomeu <
>>>>>>>>>> thibault.bridelbertomeu at gmail.com> wrote:
>>>>>>>>>>
>>>>>>>>>>> Hello Vaclav,
>>>>>>>>>>>
>>>>>>>>>>> Thank you for your quick answer !!
>>>>>>>>>>> OK so, if I need to push the group, I added :
>>>>>>>>>>>
>>>>>>>>>>> call PetscViewerHDF5PushGroup(hdf5Viewer, "/fields", ierr);
>>>>>>>>>>> CHKERRA(ierr)
>>>>>>>>>>>
>>>>>>>>>>> right after the call to PetscViewerFileSetName.
>>>>>>>>>>> The result is the same, it produces the following error :
>>>>>>>>>>>
>>>>>>>>>>> *[0]PETSC ERROR: --------------------- Error Message
>>>>>>>>>>> --------------------------------------------------------------*
>>>>>>>>>>> [0]PETSC ERROR: Unexpected data in file
>>>>>>>>>>> [0]PETSC ERROR: Global size of array in file is 105, not 25300
>>>>>>>>>>> as expected
>>>>>>>>>>> [0]PETSC ERROR: See
>>>>>>>>>>> https://www.mcs.anl.gov/petsc/documentation/faq.html for
>>>>>>>>>>> trouble shooting.
>>>>>>>>>>> [0]PETSC ERROR: Petsc Development GIT revision:
>>>>>>>>>>> v3.14.4-671-g707297fd510 GIT Date: 2021-02-24 22:50:05 +0000
>>>>>>>>>>> [0]PETSC ERROR: ../../../bin/eulerian3D on a named
>>>>>>>>>>> macbook-pro-de-thibault.home by tbridel Wed Apr 7 22:05:14 2021
>>>>>>>>>>> [0]PETSC ERROR: Configure options --with-clean=1
>>>>>>>>>>> --prefix=/Users/tbridel/Documents/1-CODES/04-PETSC/build_uns3D
>>>>>>>>>>> --with-make-np=2 --with-windows-graphics=0 --with-debugging=0
>>>>>>>>>>> --download-fblaslapack --download-mpich-shared=0 --with-x=0
>>>>>>>>>>> --with-pthread=0 --with-valgrind=0 --PETSC_ARCH=macosx_uns3D
>>>>>>>>>>> --with-fc=/usr/local/bin/mpifort --with-cc=/usr/local/bin/mpicc
>>>>>>>>>>> --with-cxx=/usr/local/bin/mpic++ --with-openmp=0 --download-hypre=yes
>>>>>>>>>>> --download-sowing=yes --download-metis=yes --download-parmetis=yes
>>>>>>>>>>> --download-triangle=yes --download-tetgen=yes --download-ctetgen=yes
>>>>>>>>>>> --download-p4est=yes --download-zlib=yes --download-c2html=yes
>>>>>>>>>>> --download-eigen=yes --download-pragmatic=yes
>>>>>>>>>>> --with-hdf5-dir=/usr/local/opt/hdf5-mpi
>>>>>>>>>>> --with-cmake-dir=/usr/local/opt/cmake
>>>>>>>>>>> --with-libtoolize=/usr/local/bin/glibtoolize
>>>>>>>>>>> --with-autoreconf=/usr/local/bin/autoreconf
>>>>>>>>>>> [0]PETSC ERROR: #1 PetscViewerHDF5ReadSizes_Private() line 114
>>>>>>>>>>> in /Users/tbridel/Documents/1-CODES/04-PETSC/src/vec/is/utils/hdf5io.c
>>>>>>>>>>> [0]PETSC ERROR: #2 PetscViewerHDF5Load() line 208 in
>>>>>>>>>>> /Users/tbridel/Documents/1-CODES/04-PETSC/src/vec/is/utils/hdf5io.c
>>>>>>>>>>> [0]PETSC ERROR: #3 VecLoad_HDF5() line 132 in
>>>>>>>>>>> /Users/tbridel/Documents/1-CODES/04-PETSC/src/vec/vec/utils/vecio.c
>>>>>>>>>>> [0]PETSC ERROR: #4 VecLoad_Default() line 257 in
>>>>>>>>>>> /Users/tbridel/Documents/1-CODES/04-PETSC/src/vec/vec/utils/vecio.c
>>>>>>>>>>> [0]PETSC ERROR: #5 VecLoad_Plex_Local() line 474 in
>>>>>>>>>>> /Users/tbridel/Documents/1-CODES/04-PETSC/src/dm/impls/plex/plex.c
>>>>>>>>>>> [0]PETSC ERROR: #6 VecLoad_Plex_HDF5_Internal() line 295 in
>>>>>>>>>>> /Users/tbridel/Documents/1-CODES/04-PETSC/src/dm/impls/plex/plexhdf5.c
>>>>>>>>>>> [0]PETSC ERROR: #7 VecLoad_Plex() line 496 in
>>>>>>>>>>> /Users/tbridel/Documents/1-CODES/04-PETSC/src/dm/impls/plex/plex.c
>>>>>>>>>>> [0]PETSC ERROR: #8 VecLoad() line 953 in
>>>>>>>>>>> /Users/tbridel/Documents/1-CODES/04-PETSC/src/vec/vec/interface/vector.c
>>>>>>>>>>> [0]PETSC ERROR: #9 User provided function() line 0 in User file
>>>>>>>>>>>
>>>>>>>>>>> Do you know where it could come from ?
>>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> I think I understand this. PETSc tries to be clever to allow you
>>>>>>>>>> to store timesteps. It gives the HDF5 array an extra dimension. Somehow
>>>>>>>>>> the Viewer has to know this. The TS does this automatically, so
>>>>>>>>>> you have an array
>>>>>>>>>>
>>>>>>>>>> GROUP "fields" {
>>>>>>>>>> DATASET "Solution" {
>>>>>>>>>> DATATYPE H5T_IEEE_F64LE
>>>>>>>>>> DATASPACE SIMPLE { ( 21, 5060, 5 ) / ( H5S_UNLIMITED,
>>>>>>>>>> 5060, 5 ) }
>>>>>>>>>> }
>>>>>>>>>> }
>>>>>>>>>>
>>>>>>>>>> which has 21 timesteps. However, when you create a brand new
>>>>>>>>>> Viewer, it does not know, and mistakenly thinks there is a single vector
>>>>>>>>>> of length 21 * 5 = 105. You can tell your reader which timestep
>>>>>>>>>> you want to extract using
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Viewer/PetscViewerHDF5SetTimestep.html
>>>>>>>>>>
>>>>>>>>>> Thanks,
>>>>>>>>>>
>>>>>>>>>> Matt
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>> As for what I am doing exactly : i am using the DMPlex with a
>>>>>>>>>>> PetscFV to solve the fluid mechanics Euler equations in 3D. The PetscFV
>>>>>>>>>>> linked to the DS of the DMPlex might be why there is a /fields before the
>>>>>>>>>>> /Solution maybe .. ?
>>>>>>>>>>>
>>>>>>>>>>> Cheers and thank you again for your help !!
>>>>>>>>>>>
>>>>>>>>>>> Thibault
>>>>>>>>>>>
>>>>>>>>>>> Le mer. 7 avr. 2021 à 10:07, Hapla Vaclav <
>>>>>>>>>>> vaclav.hapla at erdw.ethz.ch> a écrit :
>>>>>>>>>>>
>>>>>>>>>>>> Dear Thibault
>>>>>>>>>>>>
>>>>>>>>>>>> On 7 Apr 2021, at 08:18, Thibault Bridel-Bertomeu <
>>>>>>>>>>>> thibault.bridelbertomeu at gmail.com> wrote:
>>>>>>>>>>>>
>>>>>>>>>>>> Dear all,
>>>>>>>>>>>>
>>>>>>>>>>>> I have been facing a problem with VecLoad recently, even though
>>>>>>>>>>>> it seems to me I did exactly like in the examples/tutorials.
>>>>>>>>>>>>
>>>>>>>>>>>> Basically, a program writes a vector with the HDF5 writer like
>>>>>>>>>>>> this :
>>>>>>>>>>>>
>>>>>>>>>>>> call DMCreateGlobalVector(dm, sol, ierr); CHKERRA(ierr) call VecZeroEntries(X, ierr); CHKERRA(ierr) call PetscObjectSetName(X, "Solution", ierr); CHKERRA(ierr)
>>>>>>>>>>>>
>>>>>>>>>>>> < do something with X to fill it up with relevant data >
>>>>>>>>>>>>
>>>>>>>>>>>> call PetscViewerCreate(PETSC_COMM_WORLD, hdf5Viewer, ierr); CHKERRA(ierr) call PetscViewerSetType(hdf5Viewer, PETSCVIEWERHDF5, ierr); CHKERRA(ierr) call PetscViewerFileSetMode(hdf5Viewer, FILE_MODE_WRITE, ierr); CHKERRA(ierr); write(filename,'(A,I5.5,A)') "restart_", stepnum, ".h5" call PetscViewerFileSetName(hdf5Viewer, trim(filename), ierr); CHKERRA(ierr) call VecView(X, hdf5Viewer, ierr); CHKERRA(ierr) call PetscViewerDestroy(hdf5Viewer, ierr); CHKERRA(ierr)
>>>>>>>>>>>>
>>>>>>>>>>>> and the same program (but with different start-up options, say)
>>>>>>>>>>>> re-reads such a file like this :
>>>>>>>>>>>>
>>>>>>>>>>>> call DMCreateGlobalVector(dm, sol, ierr); CHKERRA(ierr) call VecZeroEntries(sol, ierr); CHKERRA(ierr) call PetscObjectSetName(sol, "Solution", ierr); CHKERRA(ierr)
>>>>>>>>>>>> call PetscViewerCreate(PETSC_COMM_WORLD, hdf5Viewer, ierr); CHKERRA(ierr) call PetscViewerSetType(hdf5Viewer, PETSCVIEWERHDF5, ierr); CHKERRA(ierr) call PetscViewerFileSetMode(hdf5Viewer, FILE_MODE_READ, ierr); CHKERRA(ierr) call PetscViewerFileSetName(hdf5Viewer, trim(restartname), ierr); CHKERRA(ierr) call VecLoad(sol, hdf5Viewer, ierr); CHKERRA(ierr) call PetscViewerDestroy(hdf5Viewer, ierr); CHKERRA(ierr)
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>> Such a dataset can be found under this link :
>>>>>>>>>>>> https://drive.google.com/file/d/1owLAx5vknNhj61_5ieAwnWOR9cmkTseL/view?usp=sharing
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>> I'm just looking at the HDF5 file. The structure is like this
>>>>>>>>>>>>
>>>>>>>>>>>> > $PETSC_DIR/$PETSC_ARCH/bin/h5dump -H restart_00020.h5
>>>>>>>>>>>>
>>>>>>>>>>>> HDF5 "restart_00020.h5" {
>>>>>>>>>>>> GROUP "/" {
>>>>>>>>>>>> GROUP "cell_fields" {
>>>>>>>>>>>> DATASET "Solution_FV solver" {
>>>>>>>>>>>> DATATYPE H5T_IEEE_F64LE
>>>>>>>>>>>> DATASPACE SIMPLE { ( 21, 3884, 5 ) / ( H5S_UNLIMITED,
>>>>>>>>>>>> 3884, 5 ) }
>>>>>>>>>>>> ATTRIBUTE "vector_field_type" {
>>>>>>>>>>>> DATATYPE H5T_STRING {
>>>>>>>>>>>> STRSIZE 7;
>>>>>>>>>>>> STRPAD H5T_STR_NULLTERM;
>>>>>>>>>>>> CSET H5T_CSET_ASCII;
>>>>>>>>>>>> CTYPE H5T_C_S1;
>>>>>>>>>>>> }
>>>>>>>>>>>> DATASPACE SCALAR
>>>>>>>>>>>> }
>>>>>>>>>>>> }
>>>>>>>>>>>> }
>>>>>>>>>>>> GROUP "fields" {
>>>>>>>>>>>> DATASET "Solution" {
>>>>>>>>>>>> DATATYPE H5T_IEEE_F64LE
>>>>>>>>>>>> DATASPACE SIMPLE { ( 21, 5060, 5 ) / ( H5S_UNLIMITED,
>>>>>>>>>>>> 5060, 5 ) }
>>>>>>>>>>>> }
>>>>>>>>>>>> }
>>>>>>>>>>>> DATASET "time" {
>>>>>>>>>>>> DATATYPE H5T_IEEE_F64LE
>>>>>>>>>>>> DATASPACE SIMPLE { ( 21, 1 ) / ( H5S_UNLIMITED, 1 ) }
>>>>>>>>>>>> }
>>>>>>>>>>>> }
>>>>>>>>>>>> }
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>> I would like the reader to read the /fields/Solution group
>>>>>>>>>>>> basically, but I am not even sure it tries to do that.
>>>>>>>>>>>> Anyhow, I got an error, saying that the size found in the file
>>>>>>>>>>>> (105) does not match the expected size (25300). If I look at the shape of
>>>>>>>>>>>> /fields/Solution it is given as (21, 5030, 5). First, it is weird, cause
>>>>>>>>>>>> the 21 seems to be 1 + current iteration number ... but anyways we find the
>>>>>>>>>>>> 5 variables and the 5030 cells. Only the reader seems to do 21 * 5 when it
>>>>>>>>>>>> should be doing 5030 * 5 ...
>>>>>>>>>>>> I tried adding 'PetscViewerHDF5PushGroup(hdf5Viewer,
>>>>>>>>>>>> "/fields/Solution", ierr)' to force it to read that group, but it does not
>>>>>>>>>>>> change anything.
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>> You definitely need to push the group, unless it's the root
>>>>>>>>>>>> group "/". There is no way the reader would guess the correct group if it's
>>>>>>>>>>>> not the root one [and I don't think it would be a good idea to implement
>>>>>>>>>>>> such searching].
>>>>>>>>>>>>
>>>>>>>>>>>> If you tried adding
>>>>>>>>>>>> PetscViewerHDF5PushGroup(hdf5Viewer, "/fields/Solution", ierr)
>>>>>>>>>>>> you likely pushed a wrong group. If the Vec name was set to
>>>>>>>>>>>> "Solution" like in your snippet [using PetscObjectSetName()], the absolute
>>>>>>>>>>>> dataset name to look up would be "/fields/Solution/Solution".
>>>>>>>>>>>>
>>>>>>>>>>>> But in your file, there's just a dataset "/fields/Solution", so
>>>>>>>>>>>> its parent group is just "/fields". So please try pushing this.
>>>>>>>>>>>>
>>>>>>>>>>>> I would gladly try to reproduce your case - perhaps the error
>>>>>>>>>>>> handling should be improved so that it would guide you into the right
>>>>>>>>>>>> direction. But it would be helpful to know exactly what you're doing - the
>>>>>>>>>>>> snippet with VecView() above should produce "/Solution" dataset but in the
>>>>>>>>>>>> file you're sending, there's "/fields/Solution".
>>>>>>>>>>>>
>>>>>>>>>>>> Note also you don't need to do VecZeroEntries() before loading
>>>>>>>>>>>> because VecLoad() fully rewrites the Vec data in memory anyway.
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>> I would appreciate it if anyone could give me pointers on this
>>>>>>>>>>>> issue ...
>>>>>>>>>>>>
>>>>>>>>>>>> Thank you very much in advance !!
>>>>>>>>>>>>
>>>>>>>>>>>> Thibault
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>> Thanks,
>>>>>>>>>>>> Vaclav
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> --
>>>>>>>>>> What most experimenters take for granted before they begin their
>>>>>>>>>> experiments is infinitely more interesting than any results to which their
>>>>>>>>>> experiments lead.
>>>>>>>>>> -- Norbert Wiener
>>>>>>>>>>
>>>>>>>>>> https://www.cse.buffalo.edu/~knepley/
>>>>>>>>>> <http://www.cse.buffalo.edu/~knepley/>
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>
>>>>>>>> --
>>>>>>>> What most experimenters take for granted before they begin their
>>>>>>>> experiments is infinitely more interesting than any results to which their
>>>>>>>> experiments lead.
>>>>>>>> -- Norbert Wiener
>>>>>>>>
>>>>>>>> https://www.cse.buffalo.edu/~knepley/
>>>>>>>> <http://www.cse.buffalo.edu/~knepley/>
>>>>>>>>
>>>>>>>
>>>>>
>>>>> --
>>>>> What most experimenters take for granted before they begin their
>>>>> experiments is infinitely more interesting than any results to which their
>>>>> experiments lead.
>>>>> -- Norbert Wiener
>>>>>
>>>>> https://www.cse.buffalo.edu/~knepley/
>>>>> <http://www.cse.buffalo.edu/~knepley/>
>>>>>
>>>>
>>
>> --
>> What most experimenters take for granted before they begin their
>> experiments is infinitely more interesting than any results to which their
>> experiments lead.
>> -- Norbert Wiener
>>
>> https://www.cse.buffalo.edu/~knepley/
>> <http://www.cse.buffalo.edu/~knepley/>
>>
>
>
--
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener
https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20210422/71c3c2c7/attachment-0001.html>
More information about the petsc-users
mailing list