<div dir="ltr">Thank you! This was *exactly* what I was looking it. It fixed the problem.<div><br></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Sun, Aug 9, 2015 at 2:31 PM, Fabian <span dir="ltr"><<a href="mailto:Fabian.Jakub@physik.uni-muenchen.de" target="_blank">Fabian.Jakub@physik.uni-muenchen.de</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
  
    
  
  <div text="#000000" bgcolor="#FFFFFF">
    If the problem is due to the rank-ordering, the following excerpt
    from the PETSc FAQ section may help:<br>
    <br>
<a href="http://www.mcs.anl.gov/petsc/documentation/faq.html#da_mpi_cart" target="_blank"><http://www.mcs.anl.gov/petsc/documentation/faq.html#da_mpi_cart></a><a name="14f13ee6bec6f959_da_mpi_cart"><br>
      <br>
      The PETSc DA object decomposes the domain differently than the
      MPI_Cart_create() command. How can one use them together?</a>
    <p> The MPI_Cart_create() first divides the mesh along the z
      direction, then the y, then the x. DMDA divides along the x, then
      y, then z. Thus, for example, rank 1 of the processes will be in a
      different part of the mesh for the two schemes. To resolve this
      you can create a new MPI communicator that you pass to
      DMDACreate() that renumbers the process ranks so that each
      physical process shares the same part of the mesh with both the
      DMDA and the MPI_Cart_create(). The code to determine the new
      numbering was provided by Rolf Kuiper. </p>
    <pre>// the numbers of processors per direction are (int) x_procs, y_procs, z_procs respectively
// (no parallelization in direction 'dir' means dir_procs = 1)

MPI_Comm NewComm;
int MPI_Rank, NewRank, x,y,z;

// get rank from MPI ordering:
MPI_Comm_rank(MPI_COMM_WORLD, &MPI_Rank);

// calculate coordinates of cpus in MPI ordering:
x = MPI_rank / (z_procs*y_procs);
y = (MPI_rank % (z_procs*y_procs)) / z_procs;
z = (MPI_rank % (z_procs*y_procs)) % z_procs;

// set new rank according to PETSc ordering:
NewRank = z*y_procs*x_procs + y*x_procs + x;

// create communicator with new ranks according to
PETSc ordering:
MPI_Comm_split(PETSC_COMM_WORLD, 1, NewRank, &NewComm);

// override the default communicator (was
MPI_COMM_WORLD as default)
PETSC_COMM_WORLD = NewComm;
    </pre><div><div class="h5">
    <br>
    <br>
    <div>On 08.08.2015 23:58, Matthew Knepley
      wrote:<br>
    </div>
    <blockquote type="cite">
      <div dir="ltr">
        <div class="gmail_extra">
          <div class="gmail_quote">On Sat, Aug 8, 2015 at 4:56 PM, Mani
            Chandra <span dir="ltr"><<a href="mailto:mc0710@gmail.com" target="_blank">mc0710@gmail.com</a>></span>
            wrote:<br>
            <blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
              <div dir="ltr">So basically one needs to correctly map
                <div><br>
                </div>
                <div>iPetsc, jPetsc -> iApplication, jApplication ?</div>
                <div><br>
                </div>
                <div>Is there is any standard way to do this? Can I get
                  petsc to automatically follow the same parallel
                  topology as the host application?</div>
              </div>
            </blockquote>
            <div><br>
            </div>
            <div>If you want to use DMDA, there is only one mapping of
              ranks, namely lexicographic. However, every structured
              grid code I have</div>
            <div>ever seen uses that mapping, perhaps with a permutation
              of the directions {x, y, z}. Thus, the user needs to map
              the directions</div>
            <div>in PETSc in the right order for the application. I am
              not sure how you would automate this seeing as it depends
              on the application.</div>
            <div><br>
            </div>
            <div>  Thanks,</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">
              <div dir="ltr">
                <div>Thanks,</div>
                <div>Mani</div>
              </div>
              <div class="gmail_extra"><br>
                <div class="gmail_quote">On Sat, Aug 8, 2015 at 3:12 PM,
                  Barry Smith <span dir="ltr"><<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>></span>
                  wrote:<br>
                  <blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><span><span><br>
                        > On Aug 8, 2015, at 3:08 PM, Mani Chandra
                        <<a href="mailto:mc0710@gmail.com" target="_blank">mc0710@gmail.com</a>>
                        wrote:<br>
                        ><br>
                        > Tried flipping the indices, I get a seg
                        fault.<br>
                        <br>
                      </span></span>  You would have to be careful in
                    exactly what you flip.  Note that the meaning of N1
                    and N2 etc would also be reversed between your code
                    and the PETSc DMDA code.<br>
                    <br>
                      I would create a tiny DMDA and put entires like 1
                    2 3 4 ... into the array so you can track where the
                    values go<br>
                    <span><font color="#888888"><br>
                          Barry<br>
                      </font></span>
                    <div>
                      <div><br>
                        ><br>
                        > On Sat, Aug 8, 2015 at 3:03 PM, Barry Smith
                        <<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>>
                        wrote:<br>
                        ><br>
                        > > On Aug 8, 2015, at 2:45 PM, Mani
                        Chandra <<a href="mailto:mc0710@gmail.com" target="_blank">mc0710@gmail.com</a>>
                        wrote:<br>
                        > ><br>
                        > > Thanks. Any suggestions for a fix?<br>
                        ><br>
                        >   Just flip the meaning of the x indices
                        and the y indices in the PETSc parts of the
                        code?<br>
                        ><br>
                        >   Also run with a very different N1 and  N2
                        (instead of equal size) to better test the code
                        coupling.<br>
                        ><br>
                        >   Barry<br>
                        ><br>
                        ><br>
                        > ><br>
                        > > Reorder the indices in
                        arrayApplication?<br>
                        > ><br>
                        > > On Sat, Aug 8, 2015 at 2:19 PM,
                        Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>>
                        wrote:<br>
                        > > On Sat, Aug 8, 2015 at 1:52 PM, Mani
                        Chandra <<a href="mailto:mc0710@gmail.com" target="_blank">mc0710@gmail.com</a>>
                        wrote:<br>
                        > > Hi,<br>
                        > ><br>
                        > > I'm having trouble interfacing petsc
                        to an application which I think is related to
                        the ordering of the nodes. Here's what I'm
                        trying to do:<br>
                        > ><br>
                        > > The application uses a structured grid
                        with a global array having dimensions N1 x N2,
                        which is then decomposed into a local array with
                        dimensions NX1 x NX2.<br>
                        > ><br>
                        > > I create a Petsc DMDA using<br>
                        > ><br>
                        > >     DMDACreate2d(MPI_COMM_WORLD,<br>
                        > >                  DM_BOUNDARY_PERIODIC,
                        DM_BOUNDARY_PERIODIC,<br>
                        > >                  DMDA_STENCIL_BOX,<br>
                        > >                  N1, N2,<br>
                        > >                  N1/NX1, N2/NX2,<br>
                        > >                  1, nghost,
                        PETSC_NULL, PETSC_NULL,<br>
                        > >                  &dmda);<br>
                        > ><br>
                        > > and then use this to create a vec:<br>
                        > ><br>
                        > >   DMCreateGlobalVector(dmda,
                        &vec);<br>
                        > ><br>
                        > > Now I copy the local contents of the
                        application array to the petsc array using the
                        following:<br>
                        > ><br>
                        > > Let i, j be the application indices
                        and iPetsc and jPetsc be petsc's indices, then:<br>
                        > ><br>
                        > > DMDAGetCorners(dmda, &iStart,
                        &jStart, &kStart,<br>
                        > >                                       
                          &iSize, &jSize, &kSize<br>
                        > >                               );<br>
                        > ><br>
                        > ><br>
                        > > double **arrayPetsc;<br>
                        > > DMDAVecGetArray(dmda, vec,
                        &arrayPetsc);<br>
                        > ><br>
                        > > for (int j=0, jPetsc=jStart; j<NX2,
                        jPetsc<jStart+jSize; j++, jPetsc++)<br>
                        > > {<br>
                        > >   for (int i=0, iPetsc=iStart;
                        i<NX1, iPetsc<iStart+iSize; i++, iPetsc++)<br>
                        > >   {<br>
                        > >      arrayPetsc[jPetsc][iPetsc] =
                        arrayApplication[j][i];<br>
                        > >   }<br>
                        > > }<br>
                        > ><br>
                        > > DMDAVecRestoreArray(dmda, vec,
                        &arrayPetsc);<br>
                        > ><br>
                        > > Now if I VecView(vec, viewer) and look
                        at the data that petsc has, it looks right when
                        run with 1 proc, but if I use 4 procs it's all
                        messed up (see attached plots).<br>
                        > ><br>
                        > > I should probably be using the AO
                        object but its not clear how. Could you help me
                        out?<br>
                        > ><br>
                        > > It looks like you have the global
                        order of processes reversed, meaning you have<br>
                        > ><br>
                        > >   1   3<br>
                        > ><br>
                        > >   0   2<br>
                        > ><br>
                        > > and it should be<br>
                        > ><br>
                        > >   2  3<br>
                        > ><br>
                        > >   0  1<br>
                        > ><br>
                        > >   Thanks,<br>
                        > ><br>
                        > >       Matt<br>
                        > ><br>
                        > > Thanks,<br>
                        > > Mani<br>
                        > > --<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<br>
                        > ><br>
                        ><br>
                        ><br>
                        <br>
                      </div>
                    </div>
                  </blockquote>
                </div>
                <br>
              </div>
            </blockquote>
          </div>
          <br>
          <br clear="all">
          <div><br>
          </div>
          -- <br>
          <div>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>
      </div>
    </blockquote>
    <br>
  </div></div></div>

</blockquote></div><br></div>