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

[Åsmund Ervik]

Matt,

I don't think there exists such a subroutine:

$ cd petsc-dev/src
$ ack "GetArrayDOF" | wc -l
51
$ ack "GetArrayDOFF90" | wc -l
0

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