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

[Åsmund Ervik]

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).

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.
>