<div dir="ltr">On Mon, Oct 7, 2013 at 7:59 AM, Åsmund Ervik <span dir="ltr"><<a href="mailto:asmund.ervik@ntnu.no" target="_blank">asmund.ervik@ntnu.no</a>></span> wrotMatt,<br><div class="gmail_extra"><div class="gmail_quote">
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<br>
I don't think there exists such a subroutine:<br>
<br>
$ cd petsc-dev/src<br>
$ ack "GetArrayDOF" | wc -l<br>
51<br>
$ ack "GetArrayDOFF90" | wc -l<br>
0<br></blockquote><div><br></div><div>That must have happened in the F90 rewrite. The DMDAVecGetArrayDOF() makes another</div><div>dimension for the dof. It would be easy to write the F90 wrapper if someone is going to use it.</div>
<div><br></div><div>   Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
Åsmund<br>
<br>
<br>
On 07. okt. 2013 14:39, Matthew Knepley wrote:<br>
> On Mon, Oct 7, 2013 at 7:30 AM, Åsmund Ervik <<a href="mailto:asmund.ervik@ntnu.no">asmund.ervik@ntnu.no</a>> wrote:<br>
><br>
>> Hi again,<br>
>><br>
>> I've started to delve into the world of DMDA. At this point I want to<br>
>> write code that goes as follows, but something tells me this may be a<br>
>> bad idea? Synopsis of the code given below the snippet.<br>
>><br>
>>   ! The "array" array is allocated like this (taken from examples):<br>
>>   ! call DMDAGetGhostCorners(SolScal,ib1,ibn,jb1,jbn,kb1,kbn,ierr)<br>
>>   ! allocate(array(0:dof-1,ib1:ibn,jb1:jbn,kb1:kbn))<br>
>><br>
>>   subroutine petsc_to_local(da,vec,array,f,dof)<br>
>>     DM, intent(in)                         :: da<br>
>>     Vec,intent(in)                         :: vec<br>
>>     PetscScalar, pointer, intent(inout)    :: array(:,:,:,:)<br>
>>     integer,intent(in)                     :: dof<br>
>>     real,intent(inout),dimension(1-stw:,1-stw:1-stw:,dof)<br>
>>     !<br>
>>     if (dof == 1) then<br>
>>       call VecGetArrayF90(da,vec,array(0,:,:,:),ierr)<br>
>>     else<br>
>>       call VecGetArrayF90(da,vec,array,ierr)<br>
>>     endif<br>
>>     call transform_petsc_us(array,f,dof)<br>
>>   end subroutine petsc_to_local<br>
>>   subroutine transform_petsc_us(array,f,dof)<br>
>>     !Note: assumed shape-array does the "coordinate transformation"<br>
>>     real, intent(in), dimension(dof,1-stw:,1-stw:,1-stw:) :: array<br>
>>     real,intent(inout),dimension(1-stw:,1-stw:1-stw:,dof) :: f<br>
>>     integer :: n<br>
>>     !<br>
>>     do n=1,dof<br>
>>       f(:,:,:,n) = array(n-1,:,:,:)<br>
>>     enddo<br>
>>   end subroutine transform_petsc_us<br>
>><br>
>> This is basically a wrapper around VecGetArrayF90 that does some things:<br>
>> it handles dof=1 or higher in the same way, it reindexes dof to start at<br>
>> 1, moves the dof from the first to the last index (these two are for<br>
>> compatibility with existing code), and perhaps most importantly it does<br>
>> a "coordinate transformation" such that the end result (in the array f<br>
>> that my code subsequently uses) is indexed from 1-stw (stw: stencil<br>
>> width) to ilmax+stw (ilmax: local index i max).<br>
>><br>
><br>
> You might want VecGetArrayDOFF90(). Its too bad F90 pointers cannot be<br>
> easily<br>
> manipulated like C or having a 1-offset would not require a copy.<br>
><br>
>     Matt<br>
><br>
><br>
>> This enables me to plug existing code straight in and manipulate the<br>
>> array f. Then I do the opposite remap and call VecRestoreArrayF90, and<br>
>> then DMLocalToGlobalBegin/End. The code between VecGet... and<br>
>> VecRestore... will run for several seconds (or longer), it includes a<br>
>> Poisson solve with KSP (I haven't figured that part out yet...)<br>
>><br>
>> Am I being stupid here, or is this a good solution? I'm guessing this<br>
>> ends up making a copy into f, so it may hurt performance somewhat, but I<br>
>> don't think that will be critical?<br>
>><br>
>> Best regards,<br>
>> Åsmund<br>
>><br>
>><br>
>> On 04. okt. 2013 22:40, Jed Brown wrote:<br>
>>> Matthew Knepley <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>> writes:<br>
>>><br>
>>>> On Fri, Oct 4, 2013 at 2:10 PM, Åsmund Ervik <<a href="mailto:asmund.ervik@ntnu.no">asmund.ervik@ntnu.no</a>><br>
>> wrote:<br>
>>>><br>
>>>>>  Is there any reason I should prefer Xdmf over Cgns? I think both use<br>
>>>>> hdf5 in the background.<br>
>>>>><br>
>>>><br>
>>>> I guess not. It was just easier for me to get Xdmf going. We have trial<br>
>>>> CGNS as well.<br>
>>><br>
>>> I don't think we have a CGNS writer, but DMPlex can read it.  Xdmf gives<br>
>>> you a little more flexibility in choosing a schema for the HDF5 file,<br>
>>> but you can put other stuff in CGNS files and Xdmf has its share of<br>
>>> shortcomings.  If you are happy with CGNS, I would say to use it and<br>
>>> report any issues you have.  We'll enhance the implementation so it can<br>
>>> do everything you need.<br>
>>><br>
>><br>
><br>
><br>
><br>
</blockquote></div><br><br clear="all"><div><br></div>-- <br>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
-- Norbert Wiener
</div></div>