<div dir="ltr"><div><div>Hi,<br><br></div>I finally found a quite simple solution using MPI_WRITE_FILE_AT, even though there may be a more efficient one. Since it is only used to write a single file with reduced dimensions, it should not be a hindrance. I finally don't use Petsc viewers neither a secondary 2D DMDA. The routine looks like this, if anyone is interested :<br><br>subroutine WriteVectorSelectK0(da,X,k0,filename,flg_open,ierr)<br> implicit none<br></div> DM :: da<br><div> Vec :: X<br> PetscErrorCode :: ierr<br> PetscViewer :: viewer<br> PetscBool :: flg_open<br> PetscScalar, pointer :: gX(:,:,:,:)<br> PetscScalar :: LocalArray(user%dof,user%mr,user%mz)<br> character(*) :: filename<br> PetscInt :: thefile<br> integer(kind=MPI_OFFSET_KIND) :: offset <br> PetscInt :: i,j<br></div><div> ! The slice to select in the Z direction<br></div><div> PetscInt :: k0<br><br></div><div> ! access the array, indexed with global indices<br></div><div> call DMDAVecGetArrayF90(da,X,gX,ierr)<br><br></div><div> ! open a file<br></div><div> call MPI_FILE_OPEN(MPI_COMM_WORLD, filename, &<br> MPI_MODE_WRONLY + MPI_MODE_CREATE, &<br> MPI_INFO_NULL, thefile, ierr)<br><br></div><div> ! only the processes which contain data on the slice k=k0 are interesting and will write to the file<br></div><div> if (k0.ge.user%phis .and. k0.le.user%phie) then<br> do j=user%zs,user%ze<br> offset = user%dof*((j-1)*user%mr+(user%rs-1))*8<br></div><div> ! Indexing of gX : one has to subtract 1 because of the C like indexing resulting from DMDAVecGetArrayF90<br></div><div> call MPI_FILE_WRITE_AT(thefile,offset,gX(:,user%rs-1,j-1,k0-1), &<br> & user%dof, MPI_DOUBLE_PRECISION, &<br> & MPI_STATUS_IGNORE, ierr)<br> end do<br> end if<br><br> call MPI_FILE_CLOSE(thefile, ierr)<br> <br></div><div> ! restore the array<br></div><div> call DMDAVecRestoreArrayF90(da,X,gX,ierr)<br><br>end subroutine WriteVectorSelectK0<br><br><br></div><div>Unfortunately, I did not find a way to avoid the loop on j, which are in principle quite inefficient. That is because when you reach the end of the local block in the first direction (user%re in my example), the next point where j is incremented by 1 is not contiguous in memory in the written file. So you have to change the offset and use a new call to MPI_FILE_WRITE_AT.<br><br></div><div>Two things that one should be careful with :<br><br></div><div>1. By default, MPI-IO does not seem to include the 4 or 8 bytes at the beginning of the file which are often added by FORTRAN (for example Petsc Viewer add 8 bytes, which include one integer about the size of the data written in the file). When you read your file outside of FORTRAN (e.g. with python), you should be careful about this difference.<br><br></div><div>2. On my machine, Petsc Viewer writes the data in big endian. However, the routine above gives me little endian (however this may be machine dependent).<br></div><div><br></div><div>Best<br><br></div><div>Timothée<br><br><br></div><div><br><br></div></div><div class="gmail_extra"><br><div class="gmail_quote">2015-08-29 17:46 GMT+09:00 <span dir="ltr"><<a href="mailto:timothee.nicolas@gmail.com" target="_blank">timothee.nicolas@gmail.com</a>></span>:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">I see. I will have a look. Regarding memory there would be no problem to put everything on process 0 I believe. If I can't figure out how to process with your routine, I will go back to the initial try with mpi_allgather.<br>
<br>
Thx<br>
<br>
Timothee<br>
<br>
Sent from my iPhone<br>
<div class="HOEnZb"><div class="h5"><br>
> On 2015/08/29, at 12:40, Barry Smith <<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>> wrote:<br>
><br>
><br>
> I wrote a routine DMDAGetRay() that pulls a 1 dimensional slice out of a 2d DMDA and puts it on process 0. It uses AOApplicationToPetsc() so is not truly scalable but perhaps you could take a look at that. Since you say "(I have lots of points in phi, so this can happen easily)" it may be ok for you to just stick the 2d slice on process 0 and then save it?<br>
><br>
> Barry<br>
><br>
> Without using AOApplicationToPetsc() or something similar yes it is in general a nightmare.<br>
><br>
><br>
>> On Aug 28, 2015, at 10:15 PM, Timothée Nicolas <<a href="mailto:timothee.nicolas@gmail.com">timothee.nicolas@gmail.com</a>> wrote:<br>
>><br>
>> Hi,<br>
>><br>
>> I have been thinking for several hours about this problem and can't find an efficient solution, however I imagine this must be possible somehow with Petsc.<br>
>><br>
>> My problem is the following :<br>
>><br>
>> I work in 3D (R,Z,Phi) which makes my data quite heavy and I don't want to save all the data of all my fields, even just once in a while. Instead, I would like to save in a binary file a slice at a given angle, say phi=0.<br>
>><br>
>> As I did not find if it's natively possible in Petsc, I considered creating a second 2D DMDA, on which I can create 2D vectors and view them with the binary viewer. So far so good. However, upon creating the 2D DMDA, naturally the distribution of processors does not correspond to the distribution of the 3D DMDA. So I was considering creating global arrays, filling them with the data of the 3D array in phi=0, then doing an MPI_allgather to give the information to all the processors, to be able to read the array and fill the 2D Petsc Vector with it. So the code would be something along the lines of :<br>
>><br>
>> PetscScalar, pointer :: gX2D(:,:,:)<br>
>> PetscScalar, pointer :: gX(:,:,:,:)<br>
>> ! LocalArray is locally filled<br>
>> ! It is transmitted to GlobalArray via MPI_Allgather<br>
>> real(8) :: LocalArray(user%dof,user%mr,user%mz)<br>
>> real(8) :: GlobalArray(user%dof,user%mr,user%mz)<br>
>><br>
>> call DMDAVecGetArrayF90(da_phi0,X2D,gX2D,ierr)<br>
>> call DMDAVecGetArrayF90(da,X,gX,ierr)<br>
>><br>
>> do k = user%phis,user%phie<br>
>> do j = user%zs,user%ze<br>
>> do i = user%rs,user%re<br>
>> do l=1,user%dof<br>
>> if (k.eq.phi_print) then<br>
>> ! Numbering obtained with DMDAGetArrayF90 differs from usual<br>
>> LocalArray(l,i,j) = gX(l-1,i-1,j-1,k-1)<br>
>> end if<br>
>> end do<br>
>> end do<br>
>> end do<br>
>> end do<br>
>><br>
>> nvals = user%dof*user%rm*user%zm<br>
>><br>
>> call MPI_AllGather(LocalArray(1,user%rs,user%zs), &<br>
>> & nvals,MPI_REAL, &<br>
>> & GlobalArray, &<br>
>> & nvals,MPI_REAL,MPI_COMM_WORLD,ierr)<br>
>><br>
>> do j = zs2D,ze2D<br>
>> do i = rs2D,re2D<br>
>> do l=1,user%dof<br>
>> gX2D(l-1,i-1,j-1) = GlobalArray(l,i,j)<br>
>> end do<br>
>> end do<br>
>> end do<br>
>><br>
>> call DMDAVecRestoreArrayF90(da_phi0,X2D,gX2D,ierr)<br>
>> call DMDAVecRestoreArrayF90(da,X,gX,ierr)<br>
>><br>
>> The problem is that MPI_allgather is not at all that simple. Exchanging array information is much more complicated that I had anticipated ! See this long post on stackoverflow :<br>
>><br>
>> <a href="http://stackoverflow.com/questions/17508647/sending-2d-arrays-in-fortran-with-mpi-gather" rel="noreferrer" target="_blank">http://stackoverflow.com/questions/17508647/sending-2d-arrays-in-fortran-with-mpi-gather</a><br>
>><br>
>> I could probably get it to work eventually, but it's pretty complicated, and I was wondering if there was not a simpler alternative I could not see. Besides, I am concerned about what could happen if the number of processors is so large that the 2D Vector gets less than 2 points per processor (I have lots of points in phi, so this can happen easily). Then Petsc would complain.<br>
>><br>
>> Does anybody have ideas ?<br>
>><br>
>> Best<br>
>><br>
>> Timothée<br>
><br>
</div></div></blockquote></div><br></div>