<div dir="ltr"><div class="gmail_default" style="font-size:small">Hi,</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">Thanks for your answers.</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">I just figured out the issues which are mainly due to the ill-conditioning of my matrix. I found the conditional number blows up when the beam is discretized into large number of elements.</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">Now, I am using the 1D bar model to solve the same problem. The good news is the solution is always accurate and stable even I discretized into 10 million elements.</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">When I run the model with both iterative solver(CG+BJACOBI/ASM) and direct solver(SUPER_LU) in parallelization, I got the following results:</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small"><table border="0" cellpadding="0" cellspacing="0" width="682" style="border-collapse:collapse;width:512pt">
<colgroup><col width="106" style="width:80pt">
<col width="64" span="9" style="width:48pt">
</colgroup><tbody><tr height="20" style="height:15pt">
<td height="20" colspan="3" width="234" style="height:15pt;width:176pt">Mesh size: 1 million unknowns</td>
<td width="64" style="width:48pt"></td>
<td width="64" style="width:48pt"></td>
<td width="64" style="width:48pt"></td>
<td width="64" style="width:48pt"></td>
<td width="64" style="width:48pt"></td>
<td width="64" style="width:48pt"></td>
<td width="64" style="width:48pt"></td>
</tr>
<tr height="20" style="height:15pt">
<td height="20" style="height:15pt">Processes</td>
<td align="right">1</td>
<td align="right">2</td>
<td align="right">4</td>
<td align="right">6</td>
<td align="right">8</td>
<td align="right">10</td>
<td align="right">12</td>
<td align="right">16</td>
<td align="right">20</td>
</tr>
<tr height="20" style="height:15pt">
<td height="20" style="height:15pt">CG+BJ</td>
<td align="right">0.36</td>
<td align="right">0.22</td>
<td align="right">0.15</td>
<td align="right">0.12</td>
<td align="right">0.11</td>
<td align="right">0.1</td>
<td align="right">0.096</td>
<td align="right">0.097</td>
<td align="right">0.099</td>
</tr>
<tr height="20" style="height:15pt">
<td height="20" style="height:15pt">CG+ASM</td>
<td align="right">0.47</td>
<td align="right">0.46</td>
<td align="right">0.267</td>
<td align="right">0.2</td>
<td align="right">0.17</td>
<td align="right">0.15</td>
<td align="right">0.145</td>
<td align="right">0.16</td>
<td align="right">0.15</td>
</tr>
<tr height="20" style="height:15pt">
<td height="20" style="height:15pt">SUPER_LU_DIST</td>
<td align="right">4.73</td>
<td align="right">5.4</td>
<td align="right">4.69</td>
<td align="right">4.58</td>
<td align="right">4.38</td>
<td align="right">4.2</td>
<td align="right">4.27</td>
<td align="right">4.28</td>
<td align="right">4.38</td>
</tr></tbody></table></div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">It seems the CG+BJ works correctly, i.e. time decreases fast with a few more processes and reach stable with many more cores.</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">However, I have some concerns about CG+ASM and SUPER_LU_DIST. The time of both two methods goes up when I use two processes compared with uniprocess. The tendency is more obvious when I use larger mesh size.</div><div class="gmail_default" style="font-size:small">I especially doubt the results of SUPER_LU_DIST in parallelism since the overall expedition is very small which is not expected.</div><div class="gmail_default" style="font-size:small">The runtime option I use for ASM pc and SUPER_LU_DIST solver is shown as below:</div><div class="gmail_default">ASM preconditioner: -pc_type asm -pc_asm_type basic<br></div><div class="gmail_default">SUPER_LU_DIST solver: -ksp_type preonly -pc_type lu -pc_factor_mat_solver_package superlu_dist</div><div class="gmail_default"><br></div><div class="gmail_default">I use same mpiexec -n np ./xxxx for all solvers.</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">Am I using them correctly? If so, is there anyway to speed up the computation further, especially for SUPER_LU_DIST?</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">Thank you very much!</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">Bests,</div><div class="gmail_default" style="font-size:small">Jinlei</div></div><div class="gmail_extra"><br><div class="gmail_quote">On Mon, Aug 1, 2016 at 2:10 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 class="">On Mon, Aug 1, 2016 at 12:52 PM, Jinlei Shen <span dir="ltr"><<a href="mailto:jshen25@jhu.edu" target="_blank">jshen25@jhu.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"><div style="font-size:small">Hi Barry,</div><div style="font-size:small"><br></div><div style="font-size:small">Thanks for your reply.</div><div style="font-size:small"><br></div><div style="font-size:small">Firstly, as you suggested, I checked my program under valgrind. The results for both sequential and parallel cases showed there are no memory errors detected.</div><div style="font-size:small"><br></div><div style="font-size:small">Second, I coded a sequential program without using PETSC to generate the global matrix of small mesh for the same problem. I then checked the matrix both from petsc(sequential and parallel) and serial code, and they are same.</div><div style="font-size:small">The way I assembled the global matrix in parallel is first distributing the nodes and elements into processes, then I loop with elements on the calling process to put the element stiffness into the global. Since the nodes and elements in cantilever beam are numbered successively, the connectivity is simple. I didn't use any partition tools to optimize mesh. It's also easy to determine the preallocation d_nnz and o_nnz since each node only connects the left and right nodes except for beginning and end, the maximum nonzeros in each row is 6. The MatSetValue process is shown as follows:</div><div><div> do iEL = idElStart, idElEnd</div><div> g_EL = (/2*iEL-1-1,2*iEL-1,2*iEL+1-1,<wbr>2*iEL+2-1/)</div><div> call MatSetValues(SG,4,g_EL,4,g_El,<wbr>SE,ADD_VALUES,ierr)</div><div> end do</div><div>where idElStart and idElEnd are the global number of first element and end element that the process owns, g_EL is the global index for DOF in element iEL, SE is the element stiffness which is same for all elements.</div><div><span style="font-family:arial,helvetica,sans-serif;font-size:12.8px">From above assembling, most of the elements are assembled within own process while there are few elements crossing two processes. </span><br></div><div><span style="font-family:arial,helvetica,sans-serif;font-size:12.8px"><br></span></div><div>The BC for my problem(cantilever under end point load) is to fix the first two DOF, so I called the <span style="font-family:arial,helvetica,sans-serif;font-size:12.8px">MatZeroRowsColumns to set the first two rows and columns into zero with diagonal equal to one, without changing the RHS.</span><br></div></div><div style="font-size:small"><br></div><div><font face="arial, helvetica, sans-serif"><span style="font-size:12.8px">Now some new issues show up :</span></font></div><div><font face="arial, helvetica, sans-serif"><span style="font-size:12.8px"><br></span></font></div><div><font face="arial, helvetica, sans-serif"><span style="font-size:12.8px">I run with </span></font><span style="font-size:12.8px">-ksp_monitor_true_<wbr>residual and -ksp_converged_reason, the monitor showed two different residues, one is the residue I can set(preconditioned, unpreconditioned, natural), the other is called true residue. </span></div><div style="font-size:small"></div><div><span style="font-size:12.8px">I initially thought the true residue is same as unpreconditioned based on definition. But it seems not true. Is it the norm of the residue (b-Ax) between computed RHS and true RHS? But, h</span><span style="font-size:12.8px">ow to understand unprecondition residue since its definition is b-Ax as well?</span></div></div></blockquote><div><br></div></span><div>It is the unpreconditioned residual. You must be misinterpreting. And we could determine exactly if you sent the output with the suggested options.</div><span class=""><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div><span style="font-size:12.8px">Can I set the true residue as my converging criteria?</span></div></div></blockquote><div><br></div></span><div>Use right preconditioning.</div><span class=""><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div><span style="font-size:12.8px">I found the accuracy of large mesh in my problem didn't necessary depend on the tolerance I set, either preconditioned or unpreconditioned, sometimes, it showed converged while the solution is not correct. But the true residue looks reflecting the true convergence very well, if the true residue is diverging, no matter what the first residue says, the results are bad!</span></div></div></blockquote><div><br></div></span><div>Yes, your preconditioner looks singular. Note that BJACOBI has an inner solver, and by default the is GMRES/ILU(0). I think</div><div>ILU(0) is really ill-conditioned for your problem.</div><span class=""><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div><span style="font-size:12.8px">For the preconditioner concerns, actually, I used BJACOBI before I sent the first email, since the JACOBI or PBJACOBI didn't even converge when the size was large.</span></div><div><span style="font-size:12.8px">But BJACOBI also didn't perform well in the paralleliztion for large mesh as posed in my last email, while it's fine for small size (below 10k elements)</span></div><div><span style="font-size:12.8px"><br></span></div><div><span style="font-size:12.8px">Yesterday, I tried the ASM with CG using the runtime option: </span>-pc_type asm -pc_asm_type basic -sub_pc_type lu (default is ilu).</div><div><span style="font-size:12.8px">For 15k elements mesh, I am now able to get the correct answer with 1-3, 6 and more processes, using either -sub_pc_type lu or ilu. </span></div></div></blockquote><div><br></div></span><div>Yes, LU works for your subdomain solver.</div><span class=""><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div><span style="font-size:12.8px">Based on all the results I have got, it shows the results varies a lot with different PC and seems ASM is better for large problem. </span></div></div></blockquote><div><br></div></span><div>Its not ASM so much as an LU subsolver that is better.</div><span class=""><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div><span style="font-size:12.8px">But what is the major factor to produce such difference between different PCs, since it's not just the issue of computational efficiency, but also the accuracy.</span></div><div><span style="font-size:12.8px">Also, I noticed for large mesh, the solution is unstable with small number of processes, for the 15k case, the solution is not correct with 4 and 5 processes, however, the solution becomes always correct with more than 6 processes. For the 50k mesh case, more processes are required to show the stability.</span></div></div></blockquote><div><br></div></span><div>Yes, partitioning is very important here. Since you do not have a good partition, you can get these wild variations.</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div><div class="h5"><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div>What do you think about this? Anything wrong? </div><div>Since the iterative solver in parallel is first computed locally(if this is correct), can it be possible that there are 'good' and 'bad' locals when dividing the global matrix, and the result from 'bad' local will contaminate the global results. But with more processes, such risk is reduced.</div><div><br></div><div>It is highly appreciated if you could give me some instruction for above questions.</div><div><br></div><div>Thank you very much.</div><div><br></div><div>Bests,</div><div>Jinlei</div><div><span style="font-size:12.8px"><br></span></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Fri, Jul 29, 2016 at 2:09 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"><br>
First run under valgrind all the cases to make sure there is not some use of uninitialized data or overwriting of data. Go to <a href="http://www.mcs.anl.gov/petsc" rel="noreferrer" target="_blank">http://www.mcs.anl.gov/petsc</a> follow the link to FAQ and search for valgrind (the web server seems to be broken at the moment).<br>
<br>
Second it is possible that your code the assembles the matrices and vectors is not correctly assembling it for either the sequential or parallel case. Hence a different number of processes could be generating a different linear system hence inconsistent results. How are you handling the parallelism? How do you know the matrix generated in parallel is identically to that sequentially?<br>
<br>
Simple preconditioners such as pbjacobi will converge slower and slower with more elements.<br>
<br>
Note that you should run with -ksp_monitor_true_residual and -ksp_converged_reason to make sure that the iterative solver is even converging. By default PETSc KSP solvers do not stop with a big error message if they do not converge so you need make sure they are always converging.<br>
<span><font color="#888888"><br>
Barry<br>
</font></span><div><div><br>
<br>
<br>
> On Jul 29, 2016, at 11:46 AM, Jinlei Shen <<a href="mailto:jshen25@jhu.edu" target="_blank">jshen25@jhu.edu</a>> wrote:<br>
><br>
> Dear PETSC developers,<br>
><br>
> Thank you for developing such a powerful tool for scientific computations.<br>
><br>
> I'm currently trying to run a simple cantilever beam FEM to test the scalability of PETSC on multi-processors. I also want to verify whether iterative solver or direct solver is more efficient for parallel large FEM problem.<br>
><br>
> Problem description, An Euler elementary cantilever beam with point load at the end along -y direction. Each node has 2 DOF (deflection and rotation)). MPIBAIJ is used with bs = 2, dnnz and onnz are determined based on the connectivity. Loop with elements in each processor to assemble the global matrix with same element stiffness matrix. The boundary condition is set using call MatZeroRowsColumns(SG,2,g_BC,<wbr>one,PETSC_NULL_OBJECT,PETSC_<wbr>NULL_OBJECT,ierr);<br>
><br>
> Based on what I have done, I find the computations work well, i.e the results are correct compared with theoretical solution, for small mesh size (small than 5000 elements) using both solvers with different numbers of processes.<br>
><br>
> However, there are several confusing issues when I increase the mesh size to 10000 and more elements with iterative solve(CG + PCBJACOBI)<br>
><br>
> 1. For 10k elements, I can get accurate solution using iterative solver with uni-processor(i.e. only one process). However, when I use 2-8 processes, it tells the linear solver converged with different iterations, but, the results are all different for different processes and erroneous. The wired thing is when I use >9 processes, the results are correct again. I am really confused by this. Could you explain me why? If my parallelization is not correct, why it works for small cases? And I check the global matrix and RHS vector and didn't see any mallocs during the process.<br>
><br>
> 2. For 30k elements, if I use one process, it says: Linear solve did not converge due to DIVERGED_INDEFINITE_PC. Does this commonly happen for large sparse matrix? If so, is there any stable solver or pc for large problem?<br>
><br>
><br>
> For parallel computing using direct solver(SUPERLU_DIST + PCLU), I can only get accuracy when the number of elements are below 5000. There must be something wrong. The way I use the superlu_dist solver is first convert MatType to AIJ, then call PCFactorSetMatSolverPackage, and change the PC to PCLU. Do I miss anything else to run SUPER_LU correctly?<br>
><br>
><br>
> I also use SUPER_LU and iterative solver(CG+PCBJACOBI) to solve the sequential version of the same problem. The results shows that iterative solver works well for <50k elements, while SUPER_LU only gets right solution below 5k elements. Can I say iterative solver is better than SUPER_LU for large problem? How can I improve the solver to copy with very large problem, such as million by million? Another thing is it's still doubtable of performance of SUPER_LU.<br>
><br>
> For the inaccuracy issue, do you think it may be due to the memory? However, there is no memory error showing during the execution.<br>
><br>
> I really appreciate someone could resolve those puzzles above for me. My goal is to replace the current SUPER_LU solver in my parallel CPFEM main program with the iterative solver using PETSC.<br>
><br>
><br>
> Please let me if you would like to see my code in detail.<br>
><br>
> Thank you very much.<br>
><br>
> Bests,<br>
> Jinlei<br>
><br>
><br>
><br>
><br>
><br>
><br>
><br>
<br>
</div></div></blockquote></div><br></div>
</blockquote></div></div></div><span class="HOEnZb"><font color="#888888"><br><br clear="all"><div><br></div>-- <br><div 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>
</font></span></div></div>
</blockquote></div><br></div>