<div dir="ltr"><br><div class="gmail_extra"><div class="gmail_quote">On Sun, Nov 3, 2013 at 7:40 AM, Gautam Bisht <span dir="ltr"><<a href="mailto:gbisht@lbl.gov" target="_blank">gbisht@lbl.gov</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">Hi Paolo,<div class="gmail_extra"><br><br><div class="gmail_quote"><div class="im">On Sun, Nov 3, 2013 at 3:53 AM, Paolo Orsini <span dir="ltr"><<a href="mailto:paolo.orsini@gmail.com" target="_blank">paolo.orsini@gmail.com</a>></span> wrote:<br>
<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"><div dir="ltr">Hi Jed,<div><br></div><div>I found out what the problem is.</div>
<div>First of all I printed all the matrices out, to be sure that MatConvert was working fine.</div><div>And that was ok: <span style="font-size:13px;color:rgb(80,0,80);font-family:arial,sans-serif">AIJ_</span><span style="font-size:13px;color:rgb(80,0,80);font-family:arial,sans-serif">mat had the same connectivity table than </span><span style="color:rgb(80,0,80);font-family:arial,sans-serif;font-size:13px">Adj_mat.</span></div>
<div>Then I compared the dual graph (stored in DualMat) computed by the Parmetis routine (MatMeshToCellGraph) and with <span style="color:rgb(80,0,80);font-family:arial,sans-serif;font-size:13px">MatMatTransposeMult. And they were quite different.</span></div>
<div> </div><div>The computation of the dual graph of the mesh is a bit more complicated than multiplying the adjacency matrix by its transpose, but not far off. With this operation, also the cells that share only one node are connected in the dual graph.</div>
<div>Instead, the minimum number of common nodes is >1 (2 in 2D probelms, 3 in 3D problems). In fact, this is an input of MatMeshToCellGraph, I should have understood this before.</div><div><br></div><div>This can be computed doing the transpose adjacency matrix (Adj_T), then doing the multiplication line by line of Adj time Adj_T, and discard the non zero entries coming from to elements that share a number of nodes less than the minimum number of common nodes imposed. I have not implement this yet, any suggestion is welcome.</div>
<div><br></div><div>I also found out that Schotch has a facility to compute a dual graph from a mesh, but not PTScotch.</div><div>Once the graph is computed, PTSchotch can load the central dual graph, and distribute it into several processors during the loading. </div>
<div>Am i right to say that PETSC is interfaced only with PTSchotch and not with Scotch?</div><div><br></div><div>To check if the PTSchotch partition works (within PFLOTRAN ), I am computing a DualMat with parmetis, saving it into a file. Then I recompile the code (using a petsc compiled with ptscotch), an load the DualMat from a file rather then forming a new one. I did a successful test when running on one processor. but I am having trouble when try on more. <br>
</div><div><br></div><div>I though the the dual graph was computed only once, even during the mpi process, instead it seems to be recomputed more than once. Not sure why.... sure i am missing something ???</div></div></blockquote>
<div><br></div></div><div>In PFLOTRAN, MatCreateMPIAdj() is called:</div><div>- Once, if unstructured grid is specified in IMPLICIT format [in unstructured_grid.F90];</div></div></div></div></blockquote><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"><div>- Twice, if unstructured grid is specified in EXPLICIT format [in unstructured_grid.F90].</div></div></div></div></blockquote><div><br></div><div>Typo on my part, EXPLICIT format is read from <b>unstructured_explicit.F90</b></div>
<div><br></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 class="gmail_extra"><div class="gmail_quote"><div>The first call creates a Adjacency matrix based on 'CONNECTIONS' data of the explicit grid; while the second call is used to create a Dual matrix.</div>
<div><br></div><div>In order to better understand the code, you could compile the code with 'ugrid_debug=1' and look at the Adj*.out/Dual*.out for the two runs:</div><div>- regression_tests/default/discretization/<a href="http://mixed_explicit.in" target="_blank">mixed_explicit.in</a></div>
<div>- regression_tests/default/discretization/<a href="http://mixed_implicit.in" target="_blank">mixed_implicit.in</a></div><span class="HOEnZb"><font color="#888888"><div><br></div><div>-Gautam.</div></font></span><div>
<div class="h5"><div><br></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"><div dir="ltr"><div> See below the original call from PFLOTRAN:</div>
<div><br></div><div>xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx</div><div>!!! this is the pass where the adjacency matrix is filled</div><div><div> num_cells_local_old = unstructured_grid%nlmax </div>
<div> allocate(local_vertices(unstructured_grid%max_nvert_per_cell* &</div><div> num_cells_local_old))</div><div> allocate(local_vertex_offset(num_cells_local_old+1))</div><div> local_vertices = 0</div>
<div> local_vertex_offset = 0</div><div> count = 0</div><div> local_vertex_offset(1) = 0</div><div> do local_id = 1, num_cells_local_old</div><div> do ivertex = 1, unstructured_grid%max_nvert_per_cell</div><div> if (unstructured_grid%cell_vertices(ivertex,local_id) == 0) exit</div>
<div> count = count + 1</div><div> ! local vertices must be zero-based for MatCreateMPIAdj; thus subtract 1</div><div> local_vertices(count) = &</div><div> unstructured_grid%cell_vertices(ivertex,local_id) - 1</div>
<div> enddo</div><div> local_vertex_offset(local_id+1) = count </div><div> enddo</div><div> </div><div> select case (unstructured_grid%grid_type)</div><div> case(TWO_DIM_GRID)</div><div> num_common_vertices = 2 ! cells must share at least this number of vertices</div>
<div> case(THREE_DIM_GRID)</div><div> num_common_vertices = 3 ! cells must share at least this number of vertices</div><div> case default</div><div> option%io_buffer = 'Grid type not recognized '</div>
<div> call printErrMsg(option)</div><div> end select</div><div><br></div><div> ! determine the global offset from 0 for cells on this rank</div><div> global_offset_old = 0</div><div> call MPI_Exscan(num_cells_local_old,global_offset_old, &</div>
<div> ONE_INTEGER_MPI,MPIU_INTEGER,MPI_SUM,option%mycomm,ierr)</div><div><br></div><div> !! the adjacency matrix is computed</div><div><div> call MatCreateMPIAdj(option%mycomm,num_cells_local_old, &<br>
</div>
<div> unstructured_grid%num_vertices_global, &</div><div> local_vertex_offset, &</div><div> local_vertices,PETSC_NULL_INTEGER,Adj_mat,ierr)</div><div>
<br></div></div><div>!! the mesh dual graph is computed<br></div><div>#if defined(PETSC_HAVE_PARMETIS)<br></div><div> call MatMeshToCellGraph(Adj_mat,num_common_vertices,Dual_mat,ierr)</div><div>#endif</div><div> call MatDestroy(Adj_mat,ierr)</div>
<div> <br></div><div> call UGridPartition(unstructured_grid,option,Dual_mat,is_new, &</div><div> num_cells_local_new)</div><div> </div><div> if (allocated(local_vertices)) deallocate(local_vertices)</div>
<div> if (allocated(local_vertex_offset)) deallocate(local_vertex_offset)</div></div><div><br></div><div>xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx</div><div><br></div><div>Paolo</div><div><br></div><div><br>
</div></div><div><div><div class="gmail_extra"><br><br><div class="gmail_quote">On Sat, Nov 2, 2013 at 3:57 AM, Jed Brown <span dir="ltr"><<a href="mailto:jedbrown@mcs.anl.gov" target="_blank">jedbrown@mcs.anl.gov</a>></span> wrote:<br>
<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"><div>Paolo Orsini <<a href="mailto:paolo.orsini@gmail.com" target="_blank">paolo.orsini@gmail.com</a>> writes:<br>
<br>
> Hi Jed,<br>
><br>
</div><div>> I believe you are right.<br>
> Multiplying the connectivity matrix by its tranpose returns an ncell x<br>
> ncell matrix with non-zero entries only for those cell sharing at least one<br>
> node. Hence the Dual Graph we need....<br>
><br>
> Following your suggestions I change the code as follows:<br>
><br>
> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!<br>
> call MatCreateMPIAdj(option%mycomm,num_cells_local_old, &<br>
> unstructured_grid%num_vertices_global, &<br>
> local_vertex_offset, &<br>
> local_vertices,PETSC_NULL_INTEGER,Adj_mat,ierr)<br>
><br>
> call MatConvert(Adj_mat, MATAIJ,MAT_INITIAL_MATRIX,AIJ_mat,ierr)<br>
><br>
> call<br>
> MatMatTransposeMult(AIJ_mat,AIJ_mat,MAT_INITIAL_MATRIX,PETSC_DEFAULT_DOUBLE_PRECISION,Dual_mat,ierr)<br>
> !! call MatMeshToCellGraph(Adj_mat,num_common_vertices,Dual_mat,ierr)<br>
><br>
> call MatDestroy(Adj_mat,ierr)<br>
> call MatDestroy(AIJ_mat,ierr)<br>
> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!<br>
><br>
> When I try to run, the MatConvert takes ages, (I suppose I can form<br>
> directly the AIJ_mat),<br>
<br>
</div>I see MatConvert_Basic does not even try to preallocate, which is<br>
certainly why it takes a long time. That's easy to fix in PETSc (at<br>
least for AIJ).<br>
<br>
PETSc devs, why doesn't it make one pass counting and another pass<br>
inserting?<br>
<div><br>
> the PTScotch partitioning seems quick, but then I get an error later<br>
> on in the code, when attempting to compute the internal<br>
> connection.... Something is wrong...<br>
<br>
</div>Always send the entire error message.<br>
<div><div><br>
> I need to double check... do you see anything wrong so far?<br>
><br>
> P<br>
><br>
><br>
><br>
><br>
><br>
><br>
> On Fri, Nov 1, 2013 at 10:18 PM, Jed Brown <<a href="mailto:jedbrown@mcs.anl.gov" target="_blank">jedbrown@mcs.anl.gov</a>> wrote:<br>
><br>
>> Please keep the discussion on at least one list so it can be archived<br>
>> and so others can comment.<br>
>><br>
>> Paolo Orsini <<a href="mailto:paolo.orsini@gmail.com" target="_blank">paolo.orsini@gmail.com</a>> writes:<br>
>><br>
>> > I really don't know how the DualGraph is going to affect the the parallel<br>
>> > computation efficiency...<br>
>> > I am trying to use this partitioning as a black box, not sure what's the<br>
>> > right thing to do..<br>
>> ><br>
>> > In the code calling the partitioning (PFLOTRAN), i can see that MPIAdj is<br>
>> > built, then the Dual graph is computed, and finally the partitioning<br>
>> takes<br>
>> > place, see below:<br>
>> ><br>
>> ><br>
>> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!<br>
>> > call MatCreateMPIAdj(option%mycomm,num_cells_local_old, &<br>
>> > unstructured_grid%num_vertices_global, &<br>
>> > local_vertex_offset, &<br>
>> > local_vertices,PETSC_NULL_INTEGER,Adj_mat,ierr)<br>
>> ><br>
>> > call MatMeshToCellGraph(Adj_mat,num_common_vertices,Dual_mat,ierr)<br>
>><br>
>> You can MatConvert Adj_mat to AIJ or assemble as AIJ to begin with, then<br>
>> use MatMatTransposeMult to create Dual_mat (which will be in AIJ<br>
>> format).<br>
>><br>
>> > .........<br>
>> ><br>
>> > call MatPartitioningCreate(option%mycomm,Part,ierr)<br>
>> > call MatPartitioningSetAdjacency(Part,Dual_mat,ierr)<br>
>> > call MatPartitioningSetFromOptions(Part,ierr)<br>
>> > call MatPartitioningApply(Part,is_new,ierr)<br>
>> > call MatPartitioningDestroy(Part,ierr)<br>
>> ><br>
>> > allocate(cell_counts(option%mycommsize))<br>
>> > call ISPartitioningCount(is_new,option%mycommsize,cell_counts,ierr)<br>
>> > num_cells_local_new = cell_counts(option%myrank+1)<br>
>> > call<br>
>> > MPI_Allreduce(num_cells_local_new,iflag,ONE_INTEGER_MPI,MPIU_INTEGER, &<br>
>> > MPI_MIN,option%mycomm,ierr)<br>
>> ><br>
>> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!<br>
>> ><br>
>> > You are telling me we can build a dual graph simply multiplying Adj_mat<br>
>> by<br>
>> > its transpose?<br>
>><br>
>> Yeah, doesn't it connect all cells that share at least one vertex? That<br>
>> is exactly the nonzero pattern of the matrix product.<br>
>><br>
>> > Why using a simple Aij when we are pre-processing a parallel process?<br>
>><br>
>> MPIAIJ implements MatMatTransposeMult, but we don't have an<br>
>> implementation for MPIAdj (without calling ParMETIS).<br>
>><br>
>> > Sorry, may be what I want to do is less black box than thought...<br>
>> ><br>
>> > And thanks again for your patience and help<br>
>> ><br>
>> > Paolo<br>
>> ><br>
>> ><br>
>> ><br>
>> ><br>
>> > On Fri, Nov 1, 2013 at 9:17 PM, Jed Brown <<a href="mailto:jedbrown@mcs.anl.gov" target="_blank">jedbrown@mcs.anl.gov</a>> wrote:<br>
>> ><br>
>> >> Paolo Orsini <<a href="mailto:paolo.orsini@gmail.com" target="_blank">paolo.orsini@gmail.com</a>> writes:<br>
>> >><br>
>> >> > Hi Jed,<br>
>> >> ><br>
>> >> > Thanks for your prompt reply.<br>
>> >> ><br>
>> >> > Is it possible to do the partitioning with PTSchotch, without<br>
>> >> > using MatMeshToCellGraph?<br>
>> >> ><br>
>> >> > Is there any other function (not using ParMetis) to compute the<br>
>> adjacency<br>
>> >> > graph needed by MatPartitioningSetAdjacency, for parallel unstructured<br>
>> >> > grids?<br>
>> >><br>
>> >> You can either compute the dual graph yourself or we can implement that<br>
>> >> function without calling ParMETIS. Actually, I believe the operation is<br>
>> >> just symbolic<br>
>> >><br>
>> >> DualGraph = E * E^T<br>
>> >><br>
>> >> where E is the element connectivity matrix. So you could build a normal<br>
>> >> (AIJ instead of MPIAdj) matrix with the connectivity, then call<br>
>> >> MatMatTransposeMult().<br>
>> >><br>
>> >> Would that work for you?<br>
>> >><br>
>><br>
</div></div></blockquote></div><br></div>
<p></p></div></div><div>
-- <br>
You received this message because you are subscribed to the Google Groups "pflotran-dev" group.<br></div>
To view this discussion on the web visit <a href="https://groups.google.com/d/msgid/pflotran-dev/CACGnotZfbdvEieXFn%2B6RYzj4o4yM_%3DK0NkHcFeK5CESL%2BovcKw%40mail.gmail.com" target="_blank">https://groups.google.com/d/msgid/pflotran-dev/CACGnotZfbdvEieXFn%2B6RYzj4o4yM_%3DK0NkHcFeK5CESL%2BovcKw%40mail.gmail.com</a>.<br>
</blockquote></div></div></div><br></div></div>
</blockquote></div><br></div></div>