[petsc-users] VecView DMDA and HDF5 - Unable to write out files properly
Abhishek G.S.
gsabhishek1ags at gmail.com
Fri Oct 15 05:11:20 CDT 2021
I finally was able to output HDF5 files to read in Paraview.
These are the following functions that make this possible,
Here I create the domain using DMDACreate3d(For 2D I just use NZ=1 - You'll
Just end up with a thickness during visualization - Accordingly set your
DELTAZ)
void write_xdmf_header(const std::string &filename){
int rank;
MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
if(rank==0){
std::size_t name_split = filename.find_last_of(".");
std::string xdmf_name = filename.substr(0,name_split) + ".xdmf";
std::ofstream fp(xdmf_name,std::ios::out);
fp << "<?xml version=\"1.0\" ?>\n";
fp << "<!DOCTYPE Xdmf SYSTEM \"Xdmf.dtd\" []>\n";
fp << "<Xdmf Version=\"2.0\"> \n\n";
fp << " <Domain>\n";
fp << " <Grid Name=\"Domain\" GridType=\"Uniform\">\n";
fp << " <Topology TopologyType=\"3DCoRectMesh\"
NumberOfElements=\""<<NZ+1<<" "<<NY+1<<" "<<NX+1<<"\"/>\n";
fp << " <Geometry GeometryType=\"ORIGIN_DXDYDZ\" >"<<"\n";
fp << " <DataItem Format=\"XML\" NumberType=\"Double\"
Dimensions=\"3\"> "<<0<<" "<<0<<" "<<0<<" </DataItem>\n";
fp << " <DataItem Format=\"XML\" NumberType=\"Double\"
Dimensions=\"3\"> "<<DELTA_X<<" "<<DELTA_Y<<" "<<DELTA_Z<<" </DataItem>\n";
fp << " </Geometry>\n\n";
fp.close();
}
}
template <typename T>
void Field<T>:: write_xdmf(const std::string &filename){
int rank;
MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
if(rank==0){
std::size_t name_split = filename.find_last_of(".");
std::string xdmf_name = filename.substr(0,name_split) + ".xdmf";
std::ofstream fp(xdmf_name,std::ios::app);
fp << " <Attribute Name=\""<<field_name<<"\" Center=\"Cell\"
AttributeType=\"Scalar\">\n";
fp << " <DataItem Format=\"HDF\" Dimensions=\"" << Nx << " "
<< Ny << " " << Nz << "\">\n";
fp << " "<<filename<<":/"<<field_name<<"\n";
fp << " </DataItem>\n";
fp << " </Attribute>\n\n";
fp.close();
}
}
void write_xdmf_footer(const std::string &filename){
int rank;
MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
if(rank==0){
std::size_t name_split = filename.find_last_of(".");
std::string xdmf_name = filename.substr(0,name_split) + ".xdmf";
std::ofstream fp(xdmf_name,std::ios::app);
fp << " </Grid>\n";
fp << " </Domain>\n";
fp << "</Xdmf>\n";
fp.close();
}
}
//For HDF file writing
template <typename T>
void Field<T>:: write_to_file(const std::string &filename){
PetscViewer viewer;
PetscViewerHDF5Open(PETSC_COMM_WORLD,filename.c_str(),FILE_MODE_APPEND,&viewer);
PetscViewerHDF5SetBaseDimension2(viewer, PETSC_FALSE);
VecView(global_vec, viewer);
PetscViewerDestroy(&viewer);
write_xdmf(filename);
}
Thanks to Matteo (most of the above are from his scripts).
On Thu, 14 Oct 2021 at 21:28, Abhishek G.S. <gsabhishek1ags at gmail.com>
wrote:
> Thanks, Matthew for the clarification/ suggestion.
> Thanks, Matteo for the scripts, I'll give this a try and get back with an
> update
>
> On Thu, 14 Oct 2021 at 19:57, Matthew Knepley <knepley at gmail.com> wrote:
>
>> On Thu, Oct 14, 2021 at 9:21 AM Matteo Semplice <
>> matteo.semplice at uninsubria.it> wrote:
>>
>>>
>>> Il 14/10/21 14:37, Matthew Knepley ha scritto:
>>>
>>> On Wed, Oct 13, 2021 at 6:30 PM Abhishek G.S. <gsabhishek1ags at gmail.com>
>>> wrote:
>>>
>>>> Hi,
>>>> I need some help with getting the file output working right.
>>>>
>>>> I am using a DMDACreate3D to initialize my DM. This is my write
>>>> function
>>>>
>>>> void write(){
>>>> PetscViewer viewer;
>>>>
>>>> PetscViewerHDF5Open(PETSC_COMM_WORLD,filename.c_str(),FILE_MODE_WRITE,&viewer);
>>>> DMDAVecRestoreArray(dm,global_vector,global_array)
>>>> VecView(global_vec, viewer);
>>>> DMDAVecGetArray(dm,global_vector,global_array);
>>>> PetscViewerDestroy(&viewer);
>>>> }
>>>>
>>>> 1) I have 2 PDE's to solve. Still, I went ahead creating a single DM
>>>> with dof=1 and creating two vectors using the DMCreateGlobalVector(). I
>>>> want to write the file out periodically. Should I perform
>>>> DMDAVecRestoreArray and DMDAVecGetArray every time is write out the
>>>> global_vector? (I know that it is just indexing the pointers and there is
>>>> no copying of values. But I am not sure)
>>>>
>>>
>>> I don't think you need the Get/RestoreArray() calls here.
>>>
>>>
>>>> 2) I am writing out to HDF5 format. I see that the vecview is supposed
>>>> to reorder the global_vector based on the DM. However, when I read the H5
>>>> files, I get an error on ViSIT and my output image becomes a 1D image
>>>> rather than a 2D/3D. What might be the reason for this ?.
>>>> Error Msg : "In domain 0, your zonal variable "avtGhostZones" has 25600
>>>> values, but it should have 160. Some values were removed to ensure VisIt
>>>> runs smoothly"
>>>> I was using a 160x160x1 DM
>>>>
>>>
>>> I do not believe we support HDF5 <--> Visit/Paraview for DMDA. The
>>> VecView() is just writing out the vector as a linear array without mesh
>>> details. For
>>> interfacing with the visualization, I think we use .vtu files. You
>>> should be able to get this effect using
>>>
>>> VecViewFromOptions(global_vec, NULL, "-vec_view");
>>>
>>> in your code, and then
>>>
>>> -vec_view vtk:sol.vtu
>>>
>>> on the command line.
>>>
>>> Hi.
>>>
>>> If you want to stick with HDF5, you can also write a XDMF file with the
>>> grid information and open that in Paraview.
>>>
>>> I am attaching some routines that I have written to do that in a solver
>>> that deals with a time dependent PDE system with 2 variables; with them I
>>> end up with a single XDMF file that Paraview can load and which contains
>>> references to all timesteps in my simulations, with each timestep being
>>> contained in an HDF5 file on its own. The idea is to call writeDomain at
>>> the beginning of the simulation, writeHDF5 for each timestep that I want to
>>> save and writeSimulationXDMF at the end. (Warning: 3D is in use, while 2D
>>> ia almost untested...)
>>>
>>> It's not the optimal solution since (1) all timesteps could be in the
>>> same HDF5 and (2) in each HDF5 i write the vectors separately and it would
>>> be better to dump the entire data in one go and interpret them as a
>>> Nx*Ny*Nz*Nvariables data from the XDMF. Nevertheless they might be a
>>> starting point for you if you wan to try this approach.
>>>
>>> You can have HDF5 put the vectors in a single array with a time
>> dimension now. Then you just alter the xdmf to point into that array. I do
>> this
>> with the unstructured code.
>>
>> Thanks,
>>
>> Matt
>>
>>> Matteo
>>>
>>>
>>>
>>>> 3) I tried using the "petsc_gen_xdmf.py" to generate the xdmf files for
>>>> use in Paraview. Here the key ["viz/geometry"] is missing. The keys present
>>>> in the output H5 file are just the two vectors I am writing and has no info
>>>> about mesh. Isn't this supposed to come automatically since the vector is
>>>> attached to the DM? How do I sort this out?
>>>>
>>>
>>> This support is for unstructured grids, DMPlex and DMForest.
>>>
>>>
>>>> 4) Can I have multiple vectors attached to the DM by
>>>> DMCreateGlobalVector() even though I created the DMDA using dof=1.
>>>>
>>>
>>> Yes.
>>>
>>> Thanks,
>>>
>>> Matt
>>>
>>>
>>>> thanks,
>>>> Abhishek
>>>>
>>>
>>>
>>> --
>>> 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/
>>> <https://eur01.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.cse.buffalo.edu%2F~knepley%2F&data=04%7C01%7Cmatteo.semplice%40uninsubria.it%7Cf27a34f639784d0231c708d98f0f7279%7C9252ed8bdffc401c86ca6237da9991fa%7C0%7C0%7C637698118789688898%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=Xb8nyv7BIl6zSBbbLFkQzluIF8s%2BOLwUQeVfjtr3q4k%3D&reserved=0>
>>>
>>> --
>>> ---
>>> Professore Associato in Analisi Numerica
>>> Dipartimento di Scienza e Alta Tecnologia
>>> Università degli Studi dell'Insubria
>>> Via Valleggio, 11 - Como
>>>
>>>
>>
>> --
>> 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/20211015/e23a51c0/attachment.html>
More information about the petsc-users
mailing list