<div dir="ltr">On Mon, Oct 7, 2013 at 7:30 AM, Åsmund Ervik <span dir="ltr"><<a href="mailto:asmund.ervik@ntnu.no" target="_blank">asmund.ervik@ntnu.no</a>></span> wrote:<br><div class="gmail_extra"><div class="gmail_quote">
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">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></blockquote><div><br></div><div>You might want VecGetArrayDOFF90(). Its too bad F90 pointers cannot be easily</div><div>manipulated like C or having a 1-offset would not require a copy.</div>
<div><br></div><div>    Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
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>> 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>
</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>