<div dir="ltr">Thanks, Matt. <div><br><div>Following your explanations, my understanding is this "If we use multiple MPI processors, the global numbering of the vertices (global domain) will be different from that with only one processor, right? ". If this is the case, will it be easy for us to check the assembled matrix from multiple processors by comparing against that of a single processor ? </div><div><br></div><div>I checked the DM details for a simple 2D mesh. How could I find (from the DM view below) the information about the global numbering of a vertex when I use 2 processors ?</div><div>_______</div><div>| | |</div><div>| __ |__ |</div><div>| | |</div><div>|___ |___|</div><div><br></div><div><div dir="ltr">DM Object: Parallel Mesh 2 MPI processes<br> type: plex<br>Parallel Mesh in 2 dimensions:<br>Supports:<br>[0] Max support size: 3<br>[0]: 2 ----> 8<br>[0]: 2 ----> 12<br>[0]: 3 ----> 8<br>[0]: 3 ----> 9<br>[0]: 3 ----> 13<br>[0]: 4 ----> 9<br>[0]: 4 ----> 14<br>[0]: 5 ----> 10<br>[0]: 5 ----> 12<br>[0]: 6 ----> 10<br>[0]: 6 ----> 11<br>[0]: 6 ----> 13<br>[0]: 7 ----> 11<br>[0]: 7 ----> 14<br>[0]: 8 ----> 0<br>[0]: 9 ----> 1<br>[0]: 10 ----> 0<br>[0]: 11 ----> 1<br>[0]: 12 ----> 0<br>[0]: 13 ----> 0<br>[0]: 13 ----> 1<br>[0]: 14 ----> 1<br>[1] Max support size: 3<br>[1]: 2 ----> 8<br>[1]: 2 ----> 12<br>[1]: 3 ----> 8<br>[1]: 3 ----> 9<br>[1]: 3 ----> 13<br>[1]: 4 ----> 9<br>[1]: 4 ----> 14<br>[1]: 5 ----> 10<br>[1]: 5 ----> 12<br>[1]: 6 ----> 10<br>[1]: 6 ----> 11<br>[1]: 6 ----> 13<br>[1]: 7 ----> 11<br>[1]: 7 ----> 14<br>[1]: 8 ----> 0<br>[1]: 9 ----> 1<br>[1]: 10 ----> 0<br>[1]: 11 ----> 1<br>[1]: 12 ----> 0<br>[1]: 13 ----> 0<br>[1]: 13 ----> 1<br>[1]: 14 ----> 1<br>Cones:<br>[0] Max cone size: 4<br>[0]: 0 <---- 8 (0)<br>[0]: 0 <---- 13 (0)<br>[0]: 0 <---- 10 (-1)<br>[0]: 0 <---- 12 (-1)<br>[0]: 1 <---- 9 (0)<br>[0]: 1 <---- 14 (0)<br>[0]: 1 <---- 11 (-1)<br>[0]: 1 <---- 13 (-1)<br>[0]: 8 <---- 2 (0)<br>[0]: 8 <---- 3 (0)<br>[0]: 9 <---- 3 (0)<br>[0]: 9 <---- 4 (0)<br>[0]: 10 <---- 5 (0)<br>[0]: 10 <---- 6 (0)<br>[0]: 11 <---- 6 (0)<br>[0]: 11 <---- 7 (0)<br>[0]: 12 <---- 2 (0)<br>[0]: 12 <---- 5 (0)<br>[0]: 13 <---- 3 (0)<br>[0]: 13 <---- 6 (0)<br>[0]: 14 <---- 4 (0)<br>[0]: 14 <---- 7 (0)<br>[1] Max cone size: 4<br>[1]: 0 <---- 8 (0)<br>[1]: 0 <---- 13 (0)<br>[1]: 0 <---- 10 (-1)<br>[1]: 0 <---- 12 (-1)<br>[1]: 1 <---- 9 (0)<br>[1]: 1 <---- 14 (0)<br>[1]: 1 <---- 11 (-1)<br>[1]: 1 <---- 13 (-1)<br>[1]: 8 <---- 2 (0)<br>[1]: 8 <---- 3 (0)<br>[1]: 9 <---- 3 (0)<br>[1]: 9 <---- 4 (0)<br>[1]: 10 <---- 5 (0)<br>[1]: 10 <---- 6 (0)<br>[1]: 11 <---- 6 (0)<br>[1]: 11 <---- 7 (0)<br>[1]: 12 <---- 2 (0)<br>[1]: 12 <---- 5 (0)<br>[1]: 13 <---- 3 (0)<br>[1]: 13 <---- 6 (0)<br>[1]: 14 <---- 4 (0)<br>[1]: 14 <---- 7 (0)<br>coordinates with 1 fields<br> field 0 with 2 components<br>Process 0:<br> ( 2) dim 2 offset 0 0. 0.<br> ( 3) dim 2 offset 2 0.5 0.<br> ( 4) dim 2 offset 4 1. 0.<br> ( 5) dim 2 offset 6 0. 0.5<br> ( 6) dim 2 offset 8 0.5 0.5<br> ( 7) dim 2 offset 10 1. 0.5<br>Process 1:<br> ( 2) dim 2 offset 0 0. 0.5<br> ( 3) dim 2 offset 2 0.5 0.5<br> ( 4) dim 2 offset 4 1. 0.5<br> ( 5) dim 2 offset 6 0. 1.<br> ( 6) dim 2 offset 8 0.5 1.<br> ( 7) dim 2 offset 10 1. 1.<br>Labels:<br>Label 'marker':<br>[0]: 2 (1)<br>[0]: 3 (1)<br>[0]: 4 (1)<br>[0]: 5 (1)<br>[0]: 7 (1)<br>[0]: 8 (1)<br>[0]: 9 (1)<br>[0]: 12 (1)<br>[0]: 14 (1)<br>[1]: 2 (1)<br>[1]: 4 (1)<br>[1]: 5 (1)<br>[1]: 6 (1)<br>[1]: 7 (1)<br>[1]: 10 (1)<br>[1]: 11 (1)<br>[1]: 12 (1)<br>[1]: 14 (1)<br>Label 'Face Sets':<br>[0]: 8 (1)<br>[0]: 9 (1)<br>[0]: 14 (2)<br>[0]: 12 (4)<br>[1]: 14 (2)<br>[1]: 10 (3)<br>[1]: 11 (3)<br>[1]: 12 (4)<br>Label 'celltype':<br>[0]: 2 (0)<br>[0]: 3 (0)<br>[0]: 4 (0)<br>[0]: 5 (0)<br>[0]: 6 (0)<br>[0]: 7 (0)<br>[0]: 8 (1)<br>[0]: 9 (1)<br>[0]: 10 (1)<br>[0]: 11 (1)<br>[0]: 12 (1)<br>[0]: 13 (1)<br>[0]: 14 (1)<br>[0]: 0 (4)<br>[0]: 1 (4)<br>[1]: 2 (0)<br>[1]: 3 (0)<br>[1]: 4 (0)<br>[1]: 5 (0)<br>[1]: 6 (0)<br>[1]: 7 (0)<br>[1]: 8 (1)<br>[1]: 9 (1)<br>[1]: 10 (1)<br>[1]: 11 (1)<br>[1]: 12 (1)<br>[1]: 13 (1)<br>[1]: 14 (1)<br>[1]: 0 (4)<br>[1]: 1 (4)<br>PetscSF Object: 2 MPI processes<br> type: basic<br> [0] Number of roots=15, leaves=5, remote ranks=1<br> [0] 5 <- (1,2)<br> [0] 6 <- (1,3)<br> [0] 7 <- (1,4)<br> [0] 10 <- (1,8)<br> [0] 11 <- (1,9)<br> [1] Number of roots=15, leaves=0, remote ranks=0<br> [0] Roots referenced by my leaves, by rank<br> [0] 1: 5 edges<br> [0] 5 <- 2<br> [0] 6 <- 3<br> [0] 7 <- 4<br> [0] 10 <- 8<br> [0] 11 <- 9<br> [1] Roots referenced by my leaves, by rank<br> MultiSF sort=rank-order<br>Vec Object: 2 MPI processes</div><div class="gmail-yj6qo gmail-ajU" style="outline:none;padding:10px 0px;width:22px;margin:2px 0px 0px"><br class="gmail-Apple-interchange-newline"></div></div><div><br></div></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Thu, May 18, 2023 at 10:02 PM Matthew Knepley <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>> wrote:<br></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 dir="ltr">On Thu, May 18, 2023 at 8:47 PM neil liu <<a href="mailto:liufield@gmail.com" target="_blank">liufield@gmail.com</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"><div dir="ltr">Thanks, Matt. I am using the following steps to build a local to global mapping. <div>Step 1) PetscSectionCreate ();</div><div> PetscSectionSetNumFields();</div><div> PetscSectionSetChart ();</div><div> //Set dof for each node</div><div> PetscSectionSetup (s);</div><div> Step 2) </div><div> PetscCall(DMGetLocalToGlobalMapping(dm, <ogm));</div> PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));<div><br></div><div>For 2D mesh from the example, <a href="https://wg-beginners.readthedocs.io/en/latest/tutorials/dm/plex/introductory_tutorial_plex.html" target="_blank">https://wg-beginners.readthedocs.io/en/latest/tutorials/dm/plex/introductory_tutorial_plex.html</a></div><div>the map from rank 0 (MPI) and rank 1(MPI) matches well with the global node ordering. </div><div><br></div><div>But When I tried a 3D gmsh, (2 hexahedra (12 nodes ) were split into 12 tetrahedra). ( 2 processors in total ) Each processor handles 9 nodes separately. When I checked the local to global mapping, It seems the mapping is not right. </div><div>For example, the mapping array is </div><div> Local Global (left: local node number; right: global node number). provided. <br></div><div> 0 0 <br> 1 3 <br> 2 5 <br> 3 1 <br> 4 6 <br> 5 8 <br> 6 2 <br> 7 9 <br> 8 11<br></div><div> But the coordinates between the local (left) and global nodes (right) are not the same, meaning they are not matching. Did I miss something? </div></div></blockquote><div><br></div><div>Plex reorders the nodes when it partitions the mesh. Add</div><div><br></div><div> DMViewFromOptions(dm, NULL, "-dm_view");</div><div><br></div><div>to your code and</div><div><br></div><div> -dm_view ::ascii_info_detail</div><div><br></div><div>to the command line. Then it will print out a mesh description that shows you the coordinates</div><div>of the vertices, and you can match them up.</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"><div dir="ltr"><div>In addition, the simple 3D gmsh file has been attached if you need. <br></div><div><b><font size="4">Gmsh file:</font></b></div><div>$MeshFormat<br>4.1 0 8<br>$EndMeshFormat<br>$PhysicalNames<br>3<br>2 27 "Excit"<br>2 28 "PEC"<br>3 29 "volume"<br>$EndPhysicalNames<br>$Entities<br>8 12 6 1<br>1 0 0 0 0 <br>2 0 0.25 0 0 <br>3 0 0.25 0.25 0 <br>4 0 0 0.25 0 <br>5 1 0 0 0 <br>6 1 0.25 0 0 <br>10 1 0.25 0.25 0 <br>14 1 0 0.25 0 <br>1 0 0 0 0 0.25 0 0 2 1 -2 <br>2 0 0.25 0 0 0.25 0.25 0 2 2 -3 <br>3 0 0 0.25 0 0.25 0.25 0 2 3 -4 <br>4 0 0 0 0 0 0.25 0 2 4 -1 <br>6 1 0 0 1 0.25 0 0 2 5 -6 <br>7 1 0.25 0 1 0.25 0.25 0 2 6 -10 <br>8 1 0 0.25 1 0.25 0.25 0 2 10 -14 <br>9 1 0 0 1 0 0.25 0 2 14 -5 <br>11 0 0 0 1 0 0 0 2 1 -5 <br>12 0 0.25 0 1 0.25 0 0 2 2 -6 <br>16 0 0.25 0.25 1 0.25 0.25 0 2 3 -10 <br>20 0 0 0.25 1 0 0.25 0 2 4 -14 <br>1 0 0 0 0 0.25 0.25 1 27 4 1 2 3 4 <br>13 0 0 0 1 0.25 0 1 28 4 1 12 -6 -11 <br>17 0 0.25 0 1 0.25 0.25 0 4 2 16 -7 -12 <br>21 0 0 0.25 1 0.25 0.25 1 28 4 3 20 -8 -16 <br>25 0 0 0 1 0 0.25 0 4 4 11 -9 -20 <br>26 1 0 0 1 0.25 0.25 1 28 4 6 7 8 9 <br>1 0 0 0 1 0.25 0.25 1 29 6 -1 26 13 17 21 25 <br>$EndEntities<br>$Nodes<br>17 12 1 12<br>0 1 0 1<br>1<br>0 0 0<br>0 2 0 1<br>2<br>0 0.25 0<br>0 3 0 1<br>3<br>0 0.25 0.25<br>0 4 0 1<br>4<br>0 0 0.25<br>0 5 0 1<br>5<br>1 0 0<br>0 6 0 1<br>6<br>1 0.25 0<br>0 10 0 1<br>7<br>1 0.25 0.25<br>0 14 0 1<br>8<br>1 0 0.25<br>1 11 0 1<br>9<br>0.5 0 0<br>1 12 0 1<br>10<br>0.5 0.25 0<br>1 16 0 1<br>11<br>0.5 0.25 0.25<br>1 20 0 1<br>12<br>0.5 0 0.25<br>2 1 0 0<br>2 13 0 0<br>2 21 0 0<br>2 26 0 0<br>3 1 0 0<br>$EndNodes<br>$Elements<br>5 24 1 24<br>2 1 2 2<br>1 1 2 4 <br>2 4 2 3 <br>2 13 2 4<br>3 9 1 2 <br>4 9 2 10 <br>5 5 9 10 <br>6 5 10 6 <br>2 21 2 4<br>7 11 3 4 <br>8 11 4 12 <br>9 7 11 12 <br>10 7 12 8 <br>2 26 2 2<br>11 5 6 8 <br>12 8 6 7 <br>3 1 4 12<br>13 1 2 4 12 <br>14 10 9 12 2 <br>15 2 9 12 1 <br>16 9 10 12 8 <br>17 6 5 8 10 <br>18 10 5 8 9 <br>19 4 2 3 11 <br>20 10 12 11 2 <br>21 2 12 11 4 <br>22 12 10 11 7 <br>23 6 8 7 10 <br>24 10 8 7 12 <br>$EndElements<br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Wed, May 17, 2023 at 7:30 PM Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>> wrote:<br></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 dir="ltr">On Wed, May 17, 2023 at 6:58 PM neil liu <<a href="mailto:liufield@gmail.com" target="_blank">liufield@gmail.com</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"><div dir="ltr">Dear Petsc developers, <div><br></div><div>I am writing my own code to calculate the FEM matrix. The following is my general framework,</div><div><br></div><div>DMPlexCreateGmsh();</div><div>MPI_Comm_rank (Petsc_comm_world, &rank);</div><div>DMPlexDistribute (.., .., &dmDist);</div><div><br></div><div>dm = dmDist;</div><div>//This can create separate dm s for different processors. (reordering.)</div><div><br></div><div>MatCreate (Petsc_comm_world, &A)</div><div>// Loop over every tetrahedral element to calculate the local matrix for each processor. Then we can get a local matrix A for each processor.</div><div><br></div><div><i>My question is : it seems we should build a global matrix B (assemble all the As for each partition) and then transfer B to KSP. KSP will do the parallelization correctly, right? </i></div></div></blockquote><div><br></div><div>I would not suggest this. The more common strategy is to assemble each element matrix directly into the</div><div>global matrix B, by mapping the cell indices directly to global indices (rather than to local indices in the matrix A). You can do this in two stages. You can create a LocalToGlobalMapping in PETSc that maps</div><div>every local index to a global index. Then you can assemble into B exactly as you would assemble into A by calling MatSetValuesLocal().</div><div><br></div><div>DMPlex handles these mappings for you automatically, but I realize that it is a large number of things to buy into.</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"><div dir="ltr"><div>If that is right, I should define a whole domain matrix B before the partitioning (MatCreate (Petsc_comm_world, &B); ), and then use localtoglobal (which petsc function should I use? Do you have any examples.) map to add A to B at the right positions (MatSetValues) ? </div><div><br></div><div>Does that make sense? </div><div><br></div><div>Thanks, </div><div><br></div><div>Xiaodong</div><div> </div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div></div>
</blockquote></div><br clear="all"><div><br></div><span>-- </span><br><div dir="ltr"><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>
</blockquote></div>
</blockquote></div><br clear="all"><div><br></div><span>-- </span><br><div dir="ltr"><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>
</blockquote></div>