[petsc-users] Writing a domain decomposition code with PETSc

[Matthew Knepley]

On Mon, Oct 7, 2013 at 7:30 AM, Åsmund Ervik <asmund.ervik at ntnu.no> 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 <knepley at gmail.com> writes:
> >
> >> On Fri, Oct 4, 2013 at 2:10 PM, Åsmund Ervik <asmund.ervik at ntnu.no>
> 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.
> >
>



-- 
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20131007/7dd2ad1d/attachment.html>