On Mon, Jul 12, 2021 at 11:40 AM Matteo Semplice <matteo.semplice@uninsubria.it> wrote:

Dear all,

    I am experimenting with hdf5+xdmf output. At https://www.xdmf.org/index.php/XDMF_Model_and_Format I read that "XDMF uses XML to store Light data and to describe the data Model. Either HDF5[3] or binary files can be used to store Heavy data. The data Format is stored redundantly in both XML and HDF5."

However, if I call DMView(dmda,hdf5viewer) and then I run h5ls or h5stat on the resulting h5 file, I see no "geometry" section in the file. How should I write the geometry to the HDF5 file?

Here below is what I have tried.

The HDF5 stuff is only implemented for DMPlex since unstructured grids need to be explicitly stored. You can usually just define the structured grid in the XML
without putting anything in the HDF5. We could write metadata so that the XML could be autogenerated, but we have not done that.

  Thanks,

     Matt

Best

    Matteo

//Setup

ierr = DMDACreate2d(PETSC_COMM_WORLD,
                              DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,
                              DMDA_STENCIL_STAR,
                              ctx.Nx,ctx.Ny, //global dim
                              PETSC_DECIDE,PETSC_DECIDE, //n proc on each dim
                              2,stWidth, //dof, stencil width
                              NULL, NULL, //n nodes per direction on each cpu
                              &(ctx.daAll));

ierr = DMDASetUniformCoordinates(ctx.daAll, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0); CHKERRQ(ierr);
ierr = DMDASetFieldName(ctx.daAll,0,"first"); CHKERRQ(ierr);
ierr = DMDASetFieldName(ctx.daAll,1,"second"); CHKERRQ(ierr);
ierr = DMDAGetLocalInfo(ctx.daAll,&ctx.daInfo); CHKERRQ(ierr);
ierr = DMCreateFieldDecomposition(ctx.daAll,NULL, NULL, &ctx.is, &ctx.daField); CHKERRQ(ierr);

//Initial data
ierr = DMCreateGlobalVector(ctx.daAll,&ctx.U0); CHKERRQ(ierr);
ierr = VecISSet(ctx.U0,ctx.is[0],1.); CHKERRQ(ierr);
ierr = VecISSet(ctx.U0,ctx.is[1],100.0); CHKERRQ(ierr);


Write mesh and Vec to hdf5:

PetscViewer viewer;
  ierr = PetscViewerHDF5Open(PETSC_COMM_WORLD,"solution.h5",FILE_MODE_WRITE,&viewer); CHKERRQ(ierr);
  ierr = DMView(ctx.daAll , viewer); CHKERRQ(ierr);  //does not output anything to solution.h5??
  ierr = VecView(ctx.U0,viewer); CHKERRQ(ierr);

ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);

Attempt to save the two fields separately::

PetscViewer viewer;
ierr = PetscViewerHDF5Open(PETSC_COMM_WORLD,"solution.h5",FILE_MODE_WRITE,&viewer); CHKERRQ(ierr);
ierr = DMView(ctx.daField[0] , viewer); CHKERRQ(ierr); //does not output anything to solution.h5??

Vec uF;
ierr = VecGetSubVector(ctx.U,ctx.is[0],&uF); CHKERRQ(ierr);
PetscObjectSetName((PetscObject) uF, "first");
ierr = VecView(uF,viewer); CHKERRQ(ierr);
ierr = VecRestoreSubVector(ctx.U,ctx.is[0],&uF); CHKERRQ(ierr);

ierr = VecGetSubVector(ctx.U,ctx.is[1],&uF); CHKERRQ(ierr);
PetscObjectSetName((PetscObject) uF, "second");
ierr = VecView(uF,viewer); CHKERRQ(ierr);
ierr = VecRestoreSubVector(ctx.U,ctx.is[1],&uF); CHKERRQ(ierr);

ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);



--
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