Dear all,     it seems it should be a fairly straighforward thing to do but I am struggling with the output of my finite difference simulation. I have tried adapting my XML ascii VTK output routines that work nicely for finite volumes (but I have issues at points where 4 CPU subdomains touch), with an HDF5 PetscViewer (for which I cannot get a correct xdmf file) and few other random attempts, but neither got me a fully satisfactory result. Rather than correcting my attempts, I am ready to start afresh, also since this project will end up with much larger meshes than those I am used to. In your experience, what is the best way to output data - associated to a DMDA with more than one scalar fields, so that each variable can be visualized independently (dof 0 to scalar field "A", dof 1 to field "B", etc) - compatible with paraview (or visit or any other free tool on linux, if it need be) - with decent scaling, i.e. Vec data should be written directly by each CPU - binary format so that files are not incredibly huge is a plus - (bonus point) in the (not so) long run this will become a method for an arbitrary  subdomain of the Cartesian grid defined by a level-set function, so a way to distinguish the points of the DMDA that are outside the physical domain would be a bonus (for example I could easily fill the Vec with NaNs for those points as long as the output/visualization can handle this). Thanks in advance! Best     Matteo