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

[Matthew Knepley]

On Mon, Oct 7, 2013 at 7:59 AM, Åsmund Ervik <asmund.ervik at ntnu.no>wrotMatt,

>
> I don't think there exists such a subroutine:
>
> $ cd petsc-dev/src
> $ ack "GetArrayDOF" | wc -l
> 51
> $ ack "GetArrayDOFF90" | wc -l
> 0
>

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.

   Matt


> Åsmund
>
>
> On 07. okt. 2013 14:39, Matthew Knepley wrote:
> > 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/d2c19db5/attachment.html>