On Mon, Jul 12, 2021 at 11:40 AM Matteo Semplice < [email protected]> 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] <https://www.hdfgroup.org/HDF5> 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 https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
