On 07. okt. 2013 15:02, Matthew Knepley wrote: > > That must have happened in the F90 rewrite. The DMDAVecGetArrayDOF() makes > another > dimension for the dof. It would be easy to write the F90 wrapper if someone > is going to use it.
My question was not really on this point. DMDAVecGetArrayF90() works with both dof=1 and dof>1, you just have to pass in different arrays, that's no biggie. The question was really about my re-mapping of coordinates from global to local. Is what I wrote a good way to do it, or are there better ways? Åsmund > > Matt > > >> Åsmund >> >> >> On 07. okt. 2013 14:39, Matthew Knepley wrote: >>> On Mon, Oct 7, 2013 at 7:30 AM, Åsmund Ervik <[email protected]> >> wrote: >>> >>>> Hi again, >>>> >>>> I've started to delve into the world of DMDA. At this point I want to >>>> write code that goes as follows, but something tells me this may be a >>>> bad idea? Synopsis of the code given below the snippet. >>>> >>>> ! The "array" array is allocated like this (taken from examples): >>>> ! call DMDAGetGhostCorners(SolScal,ib1,ibn,jb1,jbn,kb1,kbn,ierr) >>>> ! allocate(array(0:dof-1,ib1:ibn,jb1:jbn,kb1:kbn)) >>>> >>>> subroutine petsc_to_local(da,vec,array,f,dof) >>>> DM, intent(in) :: da >>>> Vec,intent(in) :: vec >>>> PetscScalar, pointer, intent(inout) :: array(:,:,:,:) >>>> integer,intent(in) :: dof >>>> real,intent(inout),dimension(1-stw:,1-stw:1-stw:,dof) >>>> ! >>>> if (dof == 1) then >>>> call VecGetArrayF90(da,vec,array(0,:,:,:),ierr) >>>> else >>>> call VecGetArrayF90(da,vec,array,ierr) >>>> endif >>>> call transform_petsc_us(array,f,dof) >>>> end subroutine petsc_to_local >>>> subroutine transform_petsc_us(array,f,dof) >>>> !Note: assumed shape-array does the "coordinate transformation" >>>> real, intent(in), dimension(dof,1-stw:,1-stw:,1-stw:) :: array >>>> real,intent(inout),dimension(1-stw:,1-stw:1-stw:,dof) :: f >>>> integer :: n >>>> ! >>>> do n=1,dof >>>> f(:,:,:,n) = array(n-1,:,:,:) >>>> enddo >>>> end subroutine transform_petsc_us >>>> >>>> This is basically a wrapper around VecGetArrayF90 that does some things: >>>> it handles dof=1 or higher in the same way, it reindexes dof to start at >>>> 1, moves the dof from the first to the last index (these two are for >>>> compatibility with existing code), and perhaps most importantly it does >>>> a "coordinate transformation" such that the end result (in the array f >>>> that my code subsequently uses) is indexed from 1-stw (stw: stencil >>>> width) to ilmax+stw (ilmax: local index i max). >>>> >>> >>> You might want VecGetArrayDOFF90(). Its too bad F90 pointers cannot be >>> easily >>> manipulated like C or having a 1-offset would not require a copy. >>> >>> Matt >>> >>> >>>> This enables me to plug existing code straight in and manipulate the >>>> array f. Then I do the opposite remap and call VecRestoreArrayF90, and >>>> then DMLocalToGlobalBegin/End. The code between VecGet... and >>>> VecRestore... will run for several seconds (or longer), it includes a >>>> Poisson solve with KSP (I haven't figured that part out yet...) >>>> >>>> Am I being stupid here, or is this a good solution? I'm guessing this >>>> ends up making a copy into f, so it may hurt performance somewhat, but I >>>> don't think that will be critical? >>>> >>>> Best regards, >>>> Åsmund >>>> >>>> >>>> On 04. okt. 2013 22:40, Jed Brown wrote: >>>>> Matthew Knepley <[email protected]> writes: >>>>> >>>>>> On Fri, Oct 4, 2013 at 2:10 PM, Åsmund Ervik <[email protected]> >>>> wrote: >>>>>> >>>>>>> Is there any reason I should prefer Xdmf over Cgns? I think both use >>>>>>> hdf5 in the background. >>>>>>> >>>>>> >>>>>> I guess not. It was just easier for me to get Xdmf going. We have >> trial >>>>>> CGNS as well. >>>>> >>>>> I don't think we have a CGNS writer, but DMPlex can read it. Xdmf >> gives >>>>> you a little more flexibility in choosing a schema for the HDF5 file, >>>>> but you can put other stuff in CGNS files and Xdmf has its share of >>>>> shortcomings. If you are happy with CGNS, I would say to use it and >>>>> report any issues you have. We'll enhance the implementation so it can >>>>> do everything you need. >>>>> >>>> >>> >>> >>> >> > > >
