[petsc-users] VecLoad from HDF5 file
Barry Smith
bsmith at petsc.dev
Thu Apr 22 17:58:01 CDT 2021
Can a Vec have multiple sequence numbers? You have a Tao loop, a TS loop, a SNES loop and a KSP loop. VecSetStep doesn't come close to covering all of that.
Note the related topic of TSHistory, SNESHistory ....
A Vec/DM could have a PetscSequenceNumber object which could have multiple named counters. Different objects could share the same one, for example the DM and the Vec TS solution might share one. A KSP solution and KSP residual vector might share one.
BTW: I don't like the word Step for this, too generic.
Barry
> On Apr 22, 2021, at 5:37 PM, Hapla Vaclav <vaclav.hapla at erdw.ethz.ch> wrote:
>
>
>
>> On 22 Apr 2021, at 23:40, Matthew Knepley <knepley at gmail.com <mailto:knepley at gmail.com>> wrote:
>>
>> On Thu, Apr 22, 2021 at 3:00 PM Hapla Vaclav <vaclav.hapla at erdw.ethz.ch <mailto: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.
>
> What about just VecSetStep? Short, apt and sounds general enough to me.
>
> Thanks,
> Vaclav
>
>>
>> Thanks,
>>
>> Matt
>>
>> Thanks,
>> Vaclav
>>
>>> On 9 Apr 2021, at 16:34, Thibault Bridel-Bertomeu <thibault.bridelbertomeu at gmail.com <mailto: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 <mailto:knepley at gmail.com>> a écrit :
>>> On Thu, Apr 8, 2021 at 9:00 AM Thibault Bridel-Bertomeu <thibault.bridelbertomeu at gmail.com <mailto: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 <mailto: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 <mailto: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 <mailto: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 <mailto: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 <mailto:knepley at gmail.com>> a écrit :
>>> On Thu, Apr 8, 2021 at 4:08 AM Thibault Bridel-Bertomeu <thibault.bridelbertomeu at gmail.com <mailto: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 <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 <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 <mailto: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 <mailto:knepley at gmail.com>> wrote:
>>>>
>>>> On Wed, Apr 7, 2021 at 4:10 PM Thibault Bridel-Bertomeu <thibault.bridelbertomeu at gmail.com <mailto: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 <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 <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 <mailto: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 <mailto: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 <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/5bdaacdf/attachment-0001.html>
More information about the petsc-users
mailing list