<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Sat, Feb 18, 2017 at 2:54 PM, Manuel Valera <span dir="ltr"><<a href="mailto:mvalera@mail.sdsu.edu" target="_blank">mvalera@mail.sdsu.edu</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">Is very easy Matthew, im porting a serial code into a parallel petsc based one and i want to do it in the most natural straightforward way, also im learning petsc at the same time. <div><br></div><div>So the original code uses staggered grid and every variable is stored in a separate array, so i created a DA for each common array of them but they differ in just +/-1 to +/-3 in each dimension.</div><div><br></div><div>So i guess i should allocate all of them in the same DA and just ignore the extra non used points for some arrays, right?</div></div></blockquote><div><br></div><div>I have seen many staggered grid implementation in PETSc, but the most successful have been those that used a single DA for everything,</div><div>and just allowed some points at the boundary to unused. I would advocate this, but before you do it I would make a detailed drawing of the</div><div>3 layers around a corner, showing all unknowns and marking them with their indices. This is invaluable when writing the equations.</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>Regards,</div><div><br></div><div>Manuel</div></div><div class="gmail_extra"><br><div class="gmail_quote">On Sat, Feb 18, 2017 at 12:49 PM, Matthew Knepley <span dir="ltr"><<a href="mailto:knepley@gmail.com" target="_blank">knepley@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"><div class="gmail_extra"><div class="gmail_quote"><span>On Sat, Feb 18, 2017 at 2:44 PM, Manuel Valera <span dir="ltr"><<a href="mailto:mvalera@mail.sdsu.edu" target="_blank">mvalera@mail.sdsu.edu</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr">thanks guys that helped a lot! <div><br></div><div>I think i got it know, i copy the code i created in case you want to suggest something or maybe use it as example...</div><div><br></div><div>I have a bonus question: Im going to operate next in several arrays at the same time, i created them using their slightly different layouts, one DA each, a peer told me different DA are not guaranteed to share the same corners, is this correct? if so, is there a way to enforce that? Thanks</div></div></blockquote><div><br></div></span><div>I do not understand. Do these DAs have the same size, meaning number of vertices in X, Y, and Z? If so, why are you</div><div>using separate arrays. You can just use 'dof' to put all your data there, and strip out subvectors using VecStride</div><div>functions, such as </div><div><br></div><div> <a href="http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecStrideScatter.html" target="_blank">http://www.mcs.anl.gov/petsc<wbr>/petsc-current/docs/manualpage<wbr>s/Vec/VecStrideScatter.html</a></div><div><br></div><div>If the DAs are different sizes, how would they ever have the same corners.</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div><div class="m_8939547535905412967h5"><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div><div>SUBROUTINE DAs(da,veczero,localv,array)</div><div><br></div><div> ! use PetscObjects, only :: ierr</div><div><br></div><div> ! Umbrella program to update and communicate the arrays in a</div><div> ! distributed fashion using the DMDA objects from PETSc.</div><div> ! Manuel Valera 1/20/17</div><div><br></div><div> ! Arguments:</div><div> ! da = DMDA array (3D) already created and setup </div><div> ! veczero = </div><div> ! globalv = </div><div> ! localv = local chunk each processor works in.</div><div> ! array = the array to be petscified. ONLY 3D ARRAYS as now. </div><div> ! arrayp = the petsc version of the array, may be not needed.</div><div> ! cta = the scatter context to translate globalv to localv in DA indices</div><div><br></div><div> Vec,intent(inout) :: veczero,localv</div><div> Vec :: natural,globalv</div><div> PetscScalar,dimension(:,:,:) :: array</div><div> PetscScalar,pointer,dimension(<wbr>:,:,:) :: arrayp</div><div> DM,intent(inout) :: da</div><div> VecScatter :: cta</div><div> PetscInt ::xind,yind,zind,xwidth,ywidt<wbr>h,zwidth</div><div><br></div><div><br></div><div><br></div><div> !Debug:</div><div> !print*,"INFO BEFORE SCATTER:" </div><div> call DAinf(da,xind,yind,zind,xwidt<wbr>h,ywidth,zwidth)</div><div> print*, SIZE(array,1),'x',SIZE(array,2<wbr>),'x',SIZE(array,3)</div><div> !</div><div><br></div><div> call DMCreateLocalVector(da,localv,<wbr>ierr)</div><div><br></div><div> !The following if-block is an attempt to scatter the arrays as vectors</div><div> !loaded from master node.</div><div><br></div><div> if (sizex/= 1)then</div><div> !PETSC Devs suggested:</div><div> !Buffer 3d array natural has diff ordering than DMDAs:</div><div> call DMDACreateNaturalVector(da,nat<wbr>ural,ierr)</div><div> call DMDAVecGetArrayF90(da,natural,<wbr>arrayp,ierr)</div><div> !**fill up veczero***:</div><div> !TODO</div><div> if(rank==0)then</div><div> arrayp = array(xind+1:xwidth,yind+1:ywi<wbr>dth,zind+1:zwidth)</div></div><div><div> endif</div><div> call DMDAVecRestoreArrayF90(da,natu<wbr>ral,arrayp,ierr)</div><div> !call DMRestoreNaturalVector(da,natu<wbr>ral,ierr) !???? not needed?</div><div><br></div><div> !Distributing array:</div><div> call VecScatterCreateToZero(natural<wbr>,cta,veczero,ierr)</div><div> call VecScatterBegin(cta,veczero,na<wbr>tural,INSERT_VALUES,SCATTER_RE<wbr>VERSE,ierr)</div><div> call VecScatterEnd(cta,veczero,natu<wbr>ral,INSERT_VALUES,SCATTER_REVE<wbr>RSE,ierr)</div><div><br></div><div><br></div><div> call DMCreateGlobalVector(da,global<wbr>v,ierr)</div><div> call DMDANaturalToGlobalBegin(da,na<wbr>tural,INSERT_VALUES,globalv,ie<wbr>rr)</div><div> call DMDANaturalToGlobalEnd(da,natu<wbr>ral,INSERT_VALUES,globalv,ierr<wbr>)</div><div><br></div><div> call DMGlobalToLocalBegin(da,global<wbr>v,INSERT_VALUES,localv,ierr)</div><div> call DMGlobalToLocalEnd(da,globalv,<wbr>INSERT_VALUES,localv,ierr)</div><div><br></div><div> else</div><div><br></div><div> print*, "This should be called in serial only."</div><div><br></div><div> call DMDAVecGetArrayF90(da,localv,a<wbr>rrayp,ierr)</div><div><br></div><div> arrayp = array</div><div><br></div><div> call DMDAVecRestoreArrayF90(da,loca<wbr>lv,arrayp,ierr)</div><div> call DMRestoreLocalVector(da,localv<wbr>,ierr)</div><div> endif</div><div><br></div><div> call VecDestroy(globalv,ierr)</div><div> call VecScatterDestroy(cta,ierr)</div><div> call VecDestroy(natural,ierr)</div><div><br></div><div><br></div><div><br></div><div> END SUBROUTINE DAs</div></div><div><br></div><div><br></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Wed, Feb 15, 2017 at 7:31 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:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><br>
You can do this in the style of<br>
<br>
VecLoad_Binary_DA()<br>
<br>
first you take your sequential vector and make it parallel in the "natural" ordering for 3d arrays<br>
<br>
DMDACreateNaturalVector(da,&na<wbr>tural);<br>
VecScatterCreateToZero(natural<wbr>,&scatter,&veczero); /* veczero is of full size on process 0 and has zero entries on all other processes*/<br>
/* fill up veczero */<br>
VecScatterBegin(scatter,vecze<wbr>ro,natural,INSERT_VALUES,SCATT<wbr>ER_REVERSE);<br>
VecScatterEnd(scatter,veczero,<wbr>natural,INSERT_VALUES,SCATTER_<wbr>REVERSE);<br>
<br>
and then move it into the PETSc DMDA parallel ordering vector with<br>
<br>
ierr = DMCreateGlobalVector(da,&xin);<wbr>CHKERRQ(ierr);<br>
ierr = DMDANaturalToGlobalBegin(da,na<wbr>tural,INSERT_VALUES,xin);CHKER<wbr>RQ(ierr);<br>
ierr = DMDANaturalToGlobalEnd(da,natu<wbr>ral,INSERT_VALUES,xin);CHKERRQ<wbr>(ierr);<br>
<div class="m_8939547535905412967m_7877043793812707701gmail-m_-8813883284620421156HOEnZb"><div class="m_8939547535905412967m_7877043793812707701gmail-m_-8813883284620421156h5"><br>
<br>
> On Feb 15, 2017, at 7:16 PM, Manuel Valera <<a href="mailto:mvalera@mail.sdsu.edu" target="_blank">mvalera@mail.sdsu.edu</a>> wrote:<br>
><br>
> Hello,<br>
><br>
> My question this time is just if there is a way to distribute a 3D array whos located at Zero rank over the processors, if possible using the DMDAs, i'm trying not to do a lot of initialization I/O in parallel.<br>
><br>
> Thanks for your time,<br>
><br>
> Manuel<br>
<br>
</div></div></blockquote></div><br></div>
</blockquote></div></div></div><br><br clear="all"><span class="HOEnZb"><font color="#888888"><span><div><br></div>-- <br><div class="m_8939547535905412967m_7877043793812707701gmail_signature">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>
</span></font></span></div></div>
</blockquote></div><br></div>
</blockquote></div><br><br clear="all"><div><br></div>-- <br><div class="gmail_signature" data-smartmail="gmail_signature">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>