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

[Åsmund Ervik]

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