Il 14/10/21 14:37, Matthew Knepley ha scritto:
On Wed, Oct 13, 2021 at 6:30 PM Abhishek G.S.
<[email protected] <mailto:[email protected]>> 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.
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
/* Output su file HDF5 (e XDMF per ParaView) */
#include <petscviewerhdf5.h>
#include "hdf5Output.h"
#include <sstream>
#include <fstream>
PetscErrorCode HDF5output::writeDomain(AppContext &ctx){
PetscErrorCode ierr;
PetscSNPrintf(hdf5DomainName,256,"%s.h5","dominio");
PetscViewer viewer;
PetscPrintf(PETSC_COMM_WORLD,"Save domain on file %s\n",hdf5DomainName);
ierr = PetscViewerHDF5Open(PETSC_COMM_WORLD,hdf5DomainName,FILE_MODE_WRITE,&viewer); CHKERRQ(ierr);
ierr = VecView(ctx.Phi,viewer); CHKERRQ(ierr);
PetscObjectSetName((PetscObject) ctx.NODETYPE, "NodeType");
ierr = VecView(ctx.NODETYPE,viewer); CHKERRQ(ierr);
ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
std::stringstream buffer;
if (ctx.rank==0){
char xdmfname[256];
PetscSNPrintf(xdmfname,256,"%s.xdmf","dominio");
std::ofstream xdmf;
xdmf.open(xdmfname,std::ios_base::out);
xdmf << "<?xml version=\"1.0\" ?>\n";
xdmf << "<Xdmf xmlns:xi=\"http://www.w3.org/2001/XInclude\"; Version=\"2.0\">\n";
xdmf << "<Domain>\n";
xdmf << " <Grid GridType=\"Collection\" CollectionType=\"Temporal\">\n";
buffer << " <Grid GridType=\"Uniform\" Name=\"domain\">\n";
buffer << " <Time Type=\"Single\" Value=\""<<0.<<"\" />\n";
switch (ctx.dim){
case 2:
buffer << " <Topology TopologyType=\"2DCoRectMesh\" Dimensions=\"" << ctx.nnx << " " << ctx.nny << "\"/>\n";
buffer << " <Geometry GeometryType=\"ORIGIN_DXDY\">\n";
buffer << " <DataItem Format=\"XML\" NumberType=\"Float\" Dimensions=\"2\">" << ctx.xmin << " " << ctx.ymin << "</DataItem>\n";
buffer << " <DataItem Format=\"XML\" NumberType=\"Float\" Dimensions=\"2\">" << ctx.dx << " " << ctx.dy << "</DataItem>\n";
buffer << " </Geometry>\n";
break;
case 3:
buffer << " <Topology TopologyType=\"3DCoRectMesh\" Dimensions=\"" << ctx.nnx << " " << ctx.nny << " " << ctx.nnz << "\"/>\n";
buffer << " <Geometry GeometryType=\"ORIGIN_DXDYDZ\">\n";
buffer << " <DataItem Format=\"XML\" NumberType=\"Float\" Dimensions=\"3\">" << ctx.xmin << " " << ctx.ymin << " " << ctx.zmin << "</DataItem>\n";
buffer << " <DataItem Format=\"XML\" NumberType=\"Float\" Dimensions=\"3\">" << ctx.dx << " " << ctx.dy << " " << ctx.dz << "</DataItem>\n";
buffer << " </Geometry>\n";
buffer << " <Attribute Name=\"Phi\" Center=\"Node\" AttributeType=\"Scalar\">\n";
buffer << " <DataItem Format=\"HDF\" Dimensions=\"" << ctx.nnx << " " << ctx.nny << " " << ctx.nnz << "\">"<<hdf5DomainName<<":/Phi</DataItem>\n";
buffer << " </Attribute>\n";
buffer << " <Attribute Name=\"NodeType\" Center=\"Node\" AttributeType=\"Scalar\">\n";
buffer << " <DataItem Format=\"HDF\" Dimensions=\"" << ctx.nnx << " " << ctx.nny << " " << ctx.nnz << "\">"<<hdf5DomainName<<":/NodeType</DataItem>\n";
buffer << " </Attribute>\n";
break;
default:
abort();
}
buffer << " </Grid>\n";
xdmf << buffer.str();
xdmf << " </Grid>\n";
xdmf << "</Domain>\n";
xdmf << "</Xdmf>\n";
xdmf.close();
}
return ierr;
}
PetscErrorCode HDF5output::writeHDF5(Vec U, PetscScalar time, bool singleXDMF){
PetscErrorCode ierr;
if (time<nextSave)
return 0;
char hdf5name[256];
PetscSNPrintf(hdf5name,256,"%s_%d.h5",basename,stepBuffer.size());
PetscViewer viewer;
PetscPrintf(PETSC_COMM_WORLD,"Save on file %s\n",hdf5name);
ierr = PetscViewerHDF5Open(PETSC_COMM_WORLD,hdf5name,FILE_MODE_WRITE,&viewer); CHKERRQ(ierr);
Vec uField;
//s
ierr = DMGetGlobalVector(ctx.daField[var::s], &uField); CHKERRQ(ierr);
ierr = VecStrideGather(U,var::s,uField,INSERT_VALUES); CHKERRQ(ierr);
PetscObjectSetName((PetscObject) uField, "S");
ierr = VecView(uField,viewer); CHKERRQ(ierr);
ierr = DMRestoreGlobalVector(ctx.daField[var::s], &uField); CHKERRQ(ierr);
//c
ierr = DMGetGlobalVector(ctx.daField[var::c], &uField); CHKERRQ(ierr);
ierr = VecStrideGather(U,var::c,uField,INSERT_VALUES); CHKERRQ(ierr);
PetscObjectSetName((PetscObject) uField, "C");
ierr = VecView(uField,viewer); CHKERRQ(ierr);
ierr = DMRestoreGlobalVector(ctx.daField[var::c], &uField); CHKERRQ(ierr);
//levelset function and NodeTypes
//ierr = VecView(ctx.Phi,viewer); CHKERRQ(ierr);
//PetscObjectSetName((PetscObject) ctx.NODETYPE, "NodeType");
//ierr = VecView(ctx.NODETYPE,viewer); CHKERRQ(ierr);
ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
std::stringstream buffer;
if (ctx.rank==0){
buffer << " <Grid GridType=\"Uniform\" Name=\"domain\">\n";
buffer << " <Time Type=\"Single\" Value=\""<<time<<"\" />\n";
switch (ctx.dim){
case 2:
buffer << " <Topology TopologyType=\"2DCoRectMesh\" Dimensions=\"" << ctx.nnx << " " << ctx.nny << "\"/>\n";
buffer << " <Geometry GeometryType=\"ORIGIN_DXDY\">\n";
buffer << " <DataItem Format=\"XML\" NumberType=\"Float\" Dimensions=\"2\">" << ctx.xmin << " " << ctx.ymin << "</DataItem>\n";
buffer << " <DataItem Format=\"XML\" NumberType=\"Float\" Dimensions=\"2\">" << ctx.dx << " " << ctx.dy << "</DataItem>\n";
buffer << " </Geometry>\n";
buffer << " <Attribute Name=\"SO2\" Center=\"Node\" AttributeType=\"Scalar\">\n";
buffer << " <DataItem Format=\"HDF\" Dimensions=\"" << ctx.nnx << " " << ctx.nny << "\">"<<hdf5name<<":/S</DataItem>\n";
buffer << " </Attribute>\n";
buffer << " <Attribute Name=\"CaCO3\" Center=\"Node\" AttributeType=\"Scalar\">\n";
buffer << " <DataItem Format=\"HDF\" Dimensions=\"" << ctx.nnx << " " << ctx.nny << "\">"<<hdf5name<<":/C</DataItem>\n";
buffer << " </Attribute>\n";
break;
case 3:
buffer << " <Topology TopologyType=\"3DCoRectMesh\" Dimensions=\"" << ctx.nnx << " " << ctx.nny << " " << ctx.nnz << "\"/>\n";
buffer << " <Geometry GeometryType=\"ORIGIN_DXDYDZ\">\n";
buffer << " <DataItem Format=\"XML\" NumberType=\"Float\" Dimensions=\"3\">" << ctx.xmin << " " << ctx.ymin << " " << ctx.zmin << "</DataItem>\n";
buffer << " <DataItem Format=\"XML\" NumberType=\"Float\" Dimensions=\"3\">" << ctx.dx << " " << ctx.dy << " " << ctx.dz << "</DataItem>\n";
buffer << " </Geometry>\n";
buffer << " <Attribute Name=\"SO2\" Center=\"Node\" AttributeType=\"Scalar\">\n";
buffer << " <DataItem Format=\"HDF\" Dimensions=\"" << ctx.nnx << " " << ctx.nny << " " << ctx.nnz << "\">"<<hdf5name<<":/S</DataItem>\n";
buffer << " </Attribute>\n";
buffer << " <Attribute Name=\"CaCO3\" Center=\"Node\" AttributeType=\"Scalar\">\n";
buffer << " <DataItem Format=\"HDF\" Dimensions=\"" << ctx.nnx << " " << ctx.nny << " " << ctx.nnz << "\">"<<hdf5name<<":/C</DataItem>\n";
buffer << " </Attribute>\n";
buffer << " <Attribute Name=\"Phi\" Center=\"Node\" AttributeType=\"Scalar\">\n";
buffer << " <DataItem Format=\"HDF\" Dimensions=\"" << ctx.nnx << " " << ctx.nny << " " << ctx.nnz << "\">"<<hdf5DomainName<<":/Phi</DataItem>\n";
buffer << " </Attribute>\n";
buffer << " <Attribute Name=\"NodeType\" Center=\"Node\" AttributeType=\"Scalar\">\n";
buffer << " <DataItem Format=\"HDF\" Dimensions=\"" << ctx.nnx << " " << ctx.nny << " " << ctx.nnz << "\">"<<hdf5DomainName<<":/NodeType</DataItem>\n";
buffer << " </Attribute>\n";
break;
default:
abort();
}
buffer << " </Grid>\n";
}
//Note: ranks>0 have empty buffers, but stepBuffer.size does increase for everybody
hdfTime tStep;
tStep.t = time;
tStep.hdfSnippet = buffer.str();
stepBuffer.push_back(tStep);
if (singleXDMF) writeLastXDMF();
nextSave += dtSave;
return 0;
}
void HDF5output::writeLastXDMF(){
if (stepBuffer.empty())
return;
if (ctx.rank==0){
char xdmfname[256];
PetscSNPrintf(xdmfname,256,"%s_%d.xdmf",basename,stepBuffer.size()-1);
std::ofstream xdmf;
xdmf.open(xdmfname,std::ios_base::out);
xdmf << "<?xml version=\"1.0\" ?>\n";
xdmf << "<Xdmf xmlns:xi=\"http://www.w3.org/2001/XInclude\"; Version=\"2.0\">\n";
xdmf << "<Domain>\n";
xdmf << " <Grid GridType=\"Collection\" CollectionType=\"Temporal\">\n";
xdmf << stepBuffer.back().hdfSnippet;
xdmf << " </Grid>\n";
xdmf << "</Domain>\n";
xdmf << "</Xdmf>\n";
xdmf.close();
}
}
void HDF5output::writeSimulationXDMF(){
if (ctx.rank==0){
char xdmfname[256];
snprintf(xdmfname,256,"%s.xdmf",basename);
std::ofstream xdmf;
xdmf.open(xdmfname,std::ios_base::out);
xdmf << "<?xml version=\"1.0\" ?>\n";
xdmf << "<Xdmf xmlns:xi=\"http://www.w3.org/2001/XInclude\"; Version=\"2.0\">\n";
xdmf << "<Domain>\n";
xdmf << " <Grid GridType=\"Collection\" CollectionType=\"Temporal\">\n";
xdmf << " <Time TimeType=\"List\">\n";
xdmf << " <DataItem Dimensions=\""<< stepBuffer.size()<<"\">\n";
for (auto it = stepBuffer.begin(); it != stepBuffer.end(); it++)
xdmf << (*it).t << " " ;
xdmf << "\n";
xdmf << " </DataItem>\n";
xdmf << " </Time>\n";
xdmf << "\n";
for (auto it = stepBuffer.begin(); it != stepBuffer.end(); it++)
xdmf << (*it).hdfSnippet << "\n";
xdmf << " </Grid>\n";
xdmf << "</Domain>\n";
xdmf << "</Xdmf>\n";
xdmf.close();
}
}
/* Output su file VTK (per ParaView) */
#ifndef __HDF5OUTPUT_H
#define __HDF5OUTPUT_H
#include "appctx.h"
#include <vector>
#include <string>
PetscErrorCode writePhi(AppContext &_ctx);
typedef struct{
PetscScalar t;
std::string hdfSnippet;
} hdfTime;
class HDF5output{
public:
HDF5output(const char * _basename, AppContext &_ctx, PetscScalar tFinal, int nSave):
ctx(_ctx),
dtSave(tFinal/nSave),
nextSave(nSave>=0?0.:INFINITY)
{strncpy(basename,_basename,250);};
//! HDF5 output of the domain (levelset and NodeTypes)
PetscErrorCode writeDomain(AppContext &ctx);
//! HDF5 output of the solution
//! with singleXDMF=true it will also write a XDMF for this timestep
PetscErrorCode writeHDF5(Vec U, PetscScalar time, bool singleXDMF=false);
//! writes XDMF for the entire simulation
void writeSimulationXDMF();
private:
AppContext &ctx;
char basename[250];
char hdf5DomainName[256];
std::vector<hdfTime> stepBuffer;
PetscScalar dtSave, nextSave;
//! writes XDMF for a single timestep (the last)
void writeLastXDMF();
};
#endif
