[petsc-users] output DMDA to hdf5 file?

[Matthew Knepley]

On Mon, Jul 12, 2021 at 11:40 AM Matteo Semplice <
matteo.semplice at 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] <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/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20210712/4d111e31/attachment.html>