<div dir="ltr"><div dir="ltr">On Thu, Feb 25, 2021 at 4:57 PM Cameron Smith <<a href="mailto:smithc11@rpi.edu">smithc11@rpi.edu</a>> wrote:<br></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">Hello,<br>
<br>
Bringing this thread back from the dead...<br>
<br>
We made progress with creation of a distributed dmplex that matches our <br>
source mesh partition and are in need of help writing values into a <br>
vector created from the dmplex object.<br>
<br>
As discussed previously, we have created a DMPlex instance using:<br>
<br>
DMPlexCreateFromCellListPetsc(...)<br>
DMGetPointSF(...)<br>
PetscSFSetGraph(...)<br>
<br>
which gives us a distribution of mesh vertices and elements in the DM <br>
object that matches the element-based partition of our unstructured mesh.<br>
<br>
We then mark mesh vertices on the geometric model boundary using <br>
DMSetLabelValue(...) and a map from our mesh vertices to dmplex points <br>
(created during dmplex definition of vtx coordinates).<br>
<br>
Following this, we create a section for vertices:<br>
<br>
> DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd); <br>
> PetscSectionCreate(PetscObjectComm((PetscObject) dm), &s);<br>
> PetscSectionSetNumFields(s, 1);<br>
> PetscSectionSetFieldComponents(s, 0, 1);<br>
> PetscSectionSetChart(s, vStart, vEnd);<br>
> for(PetscInt v = vStart; v < vEnd; ++v) {<br>
> PetscSectionSetDof(s, v, 1);<br>
> PetscSectionSetFieldDof(s, v, 0, 1);<br>
> }<br>
> PetscSectionSetUp(s);<br>
> DMSetLocalSection(dm, s);<br>
> PetscSectionDestroy(&s);<br>
> DMGetGlobalSection(dm,&s); //update the global section<br>
<br>
We then try to write values into a local Vec for the on-process vertices <br>
(roots and leaves in sf terms) and hit an ordering problem. <br>
Specifically, we make the following sequence of calls:<br>
<br>
DMGetLocalVector(dm,&bloc);<br>
VecGetArrayWrite(bloc, &bwrite);<br>
//for loop to write values to bwrite<br>
VecRestoreArrayWrite(bloc, &bwrite);<br>
DMLocalToGlobal(dm,bloc,INSERT_VALUES,b);<br>
DMRestoreLocalVector(dm,&bloc);<br></blockquote><div><br></div><div>There is an easy way to get diagnostics here. For the local vector</div><div><br></div><div> DMGetLocalSection(dm, &s);</div><div> PetscSectionGetOffset(s, v, &off);</div><div><br></div><div>will give you the offset into the array you got from VecGetArrayWrite()</div><div>for that vertex. You can get this wrapped up using DMPlexPointLocalWrite()</div><div>which is what I tend to use for this type of stuff.</div><div><br></div><div>For the global vector</div><div><br></div><div> DMGetGlobalSection(dm, &gs);</div><div> PetscSectionGetOffset(gs, v, &off);</div><div><br></div><div>will give you the offset into the portion of the global array that is stored in this process.</div><div>If you do not own the values for this vertex, the number is negative, and it is actually</div><div>-(i+1) if the index i is the valid one on the owning process.</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
Visualizing Vec 'b' in paraview, and the<br>
original mesh, tells us that the dmplex topology and geometry (the <br>
vertex coordinates) are correct but that the order we write values is <br>
wrong (not total garbage... but clearly shifted).<br></blockquote><div><br></div><div>We do not make any guarantees that global orders match local orders. However, by default</div><div>we number up global unknowns in rank order, leaving out the dofs that we not owned.</div><div><br></div><div>Does this make sense?</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
Is there anything obviously wrong in our described approach? I suspect <br>
the section creation is wrong and/or we don't understand the order of <br>
entries in the array returned by VecGetArrayWrite.<br>
<br>
Please let us know if other info is needed. We are happy to share the <br>
relevant source code.<br>
<br>
Thank-you,<br>
Cameron<br>
<br>
<br>
On 8/25/20 8:34 AM, Cameron Smith wrote:<br>
> On 8/24/20 4:57 PM, Matthew Knepley wrote:<br>
>> On Mon, Aug 24, 2020 at 4:27 PM Jed Brown <<a href="mailto:jed@jedbrown.org" target="_blank">jed@jedbrown.org</a> <br>
>> <mailto:<a href="mailto:jed@jedbrown.org" target="_blank">jed@jedbrown.org</a>>> wrote:<br>
>><br>
>> Cameron Smith <<a href="mailto:smithc11@rpi.edu" target="_blank">smithc11@rpi.edu</a> <mailto:<a href="mailto:smithc11@rpi.edu" target="_blank">smithc11@rpi.edu</a>>> writes:<br>
>><br>
>> > We made some progress with star forest creation but still have<br>
>> work to do.<br>
>> ><br>
>> > We revisited DMPlexCreateFromCellListParallelPetsc(...) and got it<br>
>> > working by sequentially partitioning the vertex coordinates across<br>
>> > processes to satisfy the 'vertexCoords' argument. Specifically,<br>
>> rank 0<br>
>> > has the coordinates for vertices with global id 0:N/P-1, rank 1 <br>
>> has<br>
>> > N/P:2*(N/P)-1, and so on (N is the total number of global<br>
>> vertices and P<br>
>> > is the number of processes).<br>
>> ><br>
>> > The consequences of the sequential partition of vertex<br>
>> coordinates in<br>
>> > subsequent solver operations is not clear. Does it make process i<br>
>> > responsible for computations and communications associated with<br>
>> global<br>
>> > vertices i*(N/P):(i+1)*(N/P)-1 ? We assumed it does and wanted<br>
>> to confirm.<br>
>><br>
>> Yeah, in the sense that the corners would be owned by the rank you<br>
>> place them on.<br>
>><br>
>> But many methods, especially high-order, perform assembly via<br>
>> non-overlapping partition of elements, in which case the<br>
>> "computations" happen where the elements are (with any required<br>
>> vertex data for the closure of those elements being sent to the rank<br>
>> handling the element).<br>
>><br>
>> Note that a typical pattern would be to create a parallel DMPlex<br>
>> with a naive distribution, then repartition/distribute it.<br>
>><br>
>><br>
>> As Jed says, CreateParallel() just makes the most naive partition of <br>
>> vertices because we have no other information. Once<br>
>> the mesh is made, you call DMPlexDistribute() again to reduce the edge <br>
>> cut.<br>
>><br>
>> Thanks,<br>
>><br>
>> Matt<br>
>><br>
> <br>
> <br>
> Thank you.<br>
> <br>
> This is being used for PIC code with low order 2d elements whose mesh is <br>
> partitioned to minimize communications during particle operations. This <br>
> partition will not be ideal for the field solve using petsc so we're <br>
> exploring alternatives that will require minimal data movement between <br>
> the two partitions. Towards that, we'll keep pursuing the SF creation.<br>
> <br>
> -Cameron<br>
> <br>
</blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><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><br></div><div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>