<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Thu, Nov 12, 2015 at 10:02 PM, David Knezevic <span dir="ltr"><<a href="mailto:david.knezevic@akselos.com" target="_blank">david.knezevic@akselos.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">On Thu, Nov 12, 2015 at 9:58 AM, David Knezevic <span dir="ltr"><<a href="mailto:david.knezevic@akselos.com" target="_blank">david.knezevic@akselos.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"><div class="gmail_extra"><div class="gmail_quote"><span>On Thu, Nov 12, 2015 at 9:50 AM, Mark Adams <span dir="ltr"><<a href="mailto:mfadams@lbl.gov" target="_blank">mfadams@lbl.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 dir="ltr"><div><div>Note, I suspect the zero pivot is coming from a coarse grid. I don't know why no-inode fixed it. N.B., this might not be deterministic.<br><br></div>If you run with -info and grep on 'GAMG' you will see the block sizes. They should be 3 on the fine grid and the 6 on the coarse grids if you set everything up correctly.</div></div></blockquote><div><br></div></span><div>OK, thanks for the info, I'll look into this further.</div><div><br></div><div>Though note that I got it to work now without needing no-inode now. The change I made was to make sure that I matched the order of function calls from the PETSc GAMG examples. The libMesh code I was using was doing things somewhat out of order, apparently.</div></div></div></div></blockquote><div><br></div><div><br></div><div>As I mentioned above, GAMG seems to be working fine for me now. Thanks for the help on this.</div><div><br></div><div>However, I wanted to ask about another 3D elasticity test case I ran in which GAMG didn't converge. The "-ksp_view -ksp_monitor" output is below. Any insights into what options I should set in order to get it to converge in this case would be appreciated. (By the way, ML worked in this case with the default options.)</div><div><br></div><div>Thanks,</div><div>David</div><div><br></div><div>--------------------------------------------</div><div><br></div><div><div> 0 KSP Residual norm 2.475892877233e-03</div><div> 1 KSP Residual norm 6.181785698923e-05</div><div> 2 KSP Residual norm 9.439050821656e-04</div></div></div></div></div></blockquote><div><br></div><div>Something very strange is happening here. CG should converge monotonically, but above it does not. What could be happening?</div><div><br></div><div> a) We could have lost orthogonality</div><div><br></div><div> This seems really unlikely in 2 iterates.</div><div><br></div><div> b) The preconditioner could be nonlinear</div><div><br></div><div> In order to test this, could you run with both GMRES and FGMRES? If the convergence is different, then something is wrong</div><div> with the preconditioner.</div><div><br></div><div> I have looked below, and it seems like this should be linear. Mark, could something be happening in GAMG?</div><div><br></div><div>Also, you should not be using CG/V-cycle. You should be using 1 iterate of FMG, -pc_mg_type full. I have just been going over</div><div>this in class. It is normal for people to solve way past discretization error, but that does not make much sense. Its not hard to</div><div>get a sense of discretization error using a manufactured solution.</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 class="gmail_extra"><div class="gmail_quote"><div><div>KSP Object: 8 MPI processes</div><div> type: cg</div><div> maximum iterations=5000</div><div> tolerances: relative=1e-12, absolute=1e-50, divergence=10000</div><div> left preconditioning</div><div> using nonzero initial guess</div><div> using PRECONDITIONED norm type for convergence test</div><div>PC Object: 8 MPI processes</div><div> type: gamg</div><div> MG: type is MULTIPLICATIVE, levels=5 cycles=v</div><div> Cycles per PCApply=1</div><div> Using Galerkin computed coarse grid matrices</div><div> GAMG specific options</div><div> Threshold for dropping small values from graph 0</div><div> AGG specific options</div><div> Symmetric graph false</div><div> Coarse grid solver -- level -------------------------------</div><div> KSP Object: (mg_coarse_) 8 MPI processes</div><div> type: gmres</div><div> GMRES: restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement</div><div> GMRES: happy breakdown tolerance 1e-30</div><div> maximum iterations=1, initial guess is zero</div><div> tolerances: relative=1e-05, absolute=1e-50, divergence=10000</div><div> left preconditioning</div><div> using NONE norm type for convergence test</div><div> PC Object: (mg_coarse_) 8 MPI processes</div><div> type: bjacobi</div><div> block Jacobi: number of blocks = 8</div><div> Local solve is same for all blocks, in the following KSP and PC objects:</div><div> KSP Object: (mg_coarse_sub_) 1 MPI processes</div><div> type: preonly</div><div> maximum iterations=1, initial guess is zero</div><div> tolerances: relative=1e-05, absolute=1e-50, divergence=10000</div><div> left preconditioning</div><div> using NONE norm type for convergence test</div><div> PC Object: (mg_coarse_sub_) 1 MPI processes</div><div> type: lu</div><div> LU: out-of-place factorization</div><div> tolerance for zero pivot 2.22045e-14</div><div> using diagonal shift on blocks to prevent zero pivot [INBLOCKS]</div><div> matrix ordering: nd</div><div> factor fill ratio given 5, needed 1</div><div> Factored matrix follows:</div><div> Mat Object: 1 MPI processes</div><div> type: seqaij</div><div> rows=6, cols=6, bs=6</div><div> package used to perform factorization: petsc</div><div> total: nonzeros=36, allocated nonzeros=36</div><div> total number of mallocs used during MatSetValues calls =0</div><div> using I-node routines: found 2 nodes, limit used is 5</div><div> linear system matrix = precond matrix:</div><div> Mat Object: 1 MPI processes</div><div> type: seqaij</div><div> rows=6, cols=6, bs=6</div><div> total: nonzeros=36, allocated nonzeros=36</div><div> total number of mallocs used during MatSetValues calls =0</div><div> using I-node routines: found 2 nodes, limit used is 5</div><div> linear system matrix = precond matrix:</div><div> Mat Object: 8 MPI processes</div><div> type: mpiaij</div><div> rows=6, cols=6, bs=6</div><div> total: nonzeros=36, allocated nonzeros=36</div><div> total number of mallocs used during MatSetValues calls =0</div><div> using I-node (on process 0) routines: found 2 nodes, limit used is 5</div><div> Down solver (pre-smoother) on level 1 -------------------------------</div><div> KSP Object: (mg_levels_1_) 8 MPI processes</div><div> type: chebyshev</div><div> Chebyshev: eigenvalue estimates: min = 0.146922, max = 1.61615</div><div> Chebyshev: eigenvalues estimated using gmres with translations [0 0.1; 0 1.1]</div><div> KSP Object: (mg_levels_1_esteig_) 8 MPI processes</div><div> type: gmres</div><div> GMRES: restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement</div><div> GMRES: happy breakdown tolerance 1e-30</div><div> maximum iterations=10, initial guess is zero</div><div> tolerances: relative=1e-05, absolute=1e-50, divergence=10000</div><div> left preconditioning</div><div> using NONE norm type for convergence test</div><div> maximum iterations=2</div><div> tolerances: relative=1e-05, absolute=1e-50, divergence=10000</div><div> left preconditioning</div><div> using nonzero initial guess</div><div> using NONE norm type for convergence test</div><div> PC Object: (mg_levels_1_) 8 MPI processes</div><div> type: sor</div><div> SOR: type = local_symmetric, iterations = 1, local iterations = 1, omega = 1</div><div> linear system matrix = precond matrix:</div><div> Mat Object: 8 MPI processes</div><div> type: mpiaij</div><div> rows=162, cols=162, bs=6</div><div> total: nonzeros=26244, allocated nonzeros=26244</div><div> total number of mallocs used during MatSetValues calls =0</div><div> using I-node (on process 0) routines: found 5 nodes, limit used is 5</div><div> Up solver (post-smoother) same as down solver (pre-smoother)</div><div> Down solver (pre-smoother) on level 2 -------------------------------</div><div> KSP Object: (mg_levels_2_) 8 MPI processes</div><div> type: chebyshev</div><div> Chebyshev: eigenvalue estimates: min = 0.143941, max = 1.58335</div><div> Chebyshev: eigenvalues estimated using gmres with translations [0 0.1; 0 1.1]</div><div> KSP Object: (mg_levels_2_esteig_) 8 MPI processes</div><div> type: gmres</div><div> GMRES: restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement</div><div> GMRES: happy breakdown tolerance 1e-30</div><div> maximum iterations=10, initial guess is zero</div><div> tolerances: relative=1e-05, absolute=1e-50, divergence=10000</div><div> left preconditioning</div><div> using NONE norm type for convergence test</div><div> maximum iterations=2</div><div> tolerances: relative=1e-05, absolute=1e-50, divergence=10000</div><div> left preconditioning</div><div> using nonzero initial guess</div><div> using NONE norm type for convergence test</div><div> PC Object: (mg_levels_2_) 8 MPI processes</div><div> type: sor</div><div> SOR: type = local_symmetric, iterations = 1, local iterations = 1, omega = 1</div><div> linear system matrix = precond matrix:</div><div> Mat Object: 8 MPI processes</div><div> type: mpiaij</div><div> rows=4668, cols=4668, bs=6</div><div> total: nonzeros=3.00254e+06, allocated nonzeros=3.00254e+06</div><div> total number of mallocs used during MatSetValues calls =0</div><div> using I-node (on process 0) routines: found 158 nodes, limit used is 5</div><div> Up solver (post-smoother) same as down solver (pre-smoother)</div><div> Down solver (pre-smoother) on level 3 -------------------------------</div><div> KSP Object: (mg_levels_3_) 8 MPI processes</div><div> type: chebyshev</div><div> Chebyshev: eigenvalue estimates: min = 0.128239, max = 1.41063</div><div> Chebyshev: eigenvalues estimated using gmres with translations [0 0.1; 0 1.1]</div><div> KSP Object: (mg_levels_3_esteig_) 8 MPI processes</div><div> type: gmres</div><div> GMRES: restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement</div><div> GMRES: happy breakdown tolerance 1e-30</div><div> maximum iterations=10, initial guess is zero</div><div> tolerances: relative=1e-05, absolute=1e-50, divergence=10000</div><div> left preconditioning</div><div> using NONE norm type for convergence test</div><div> maximum iterations=2</div><div> tolerances: relative=1e-05, absolute=1e-50, divergence=10000</div><div> left preconditioning</div><div> using nonzero initial guess</div><div> using NONE norm type for convergence test</div><div> PC Object: (mg_levels_3_) 8 MPI processes</div><div> type: sor</div><div> SOR: type = local_symmetric, iterations = 1, local iterations = 1, omega = 1</div><div> linear system matrix = precond matrix:</div><div> Mat Object: 8 MPI processes</div><div> type: mpiaij</div><div> rows=66390, cols=66390, bs=6</div><div> total: nonzeros=1.65982e+07, allocated nonzeros=1.65982e+07</div><div> total number of mallocs used during MatSetValues calls =0</div><div> using I-node (on process 0) routines: found 2668 nodes, limit used is 5</div><div> Up solver (post-smoother) same as down solver (pre-smoother)</div><div> Down solver (pre-smoother) on level 4 -------------------------------</div><div> KSP Object: (mg_levels_4_) 8 MPI processes</div><div> type: chebyshev</div><div> Chebyshev: eigenvalue estimates: min = 0.1, max = 1.1</div><div> Chebyshev: eigenvalues estimated using gmres with translations [0 0.1; 0 1.1]</div><div> KSP Object: (mg_levels_4_esteig_) 8 MPI processes</div><div> type: gmres</div><div> GMRES: restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement</div><div> GMRES: happy breakdown tolerance 1e-30</div><div> maximum iterations=10, initial guess is zero</div><div> tolerances: relative=1e-05, absolute=1e-50, divergence=10000</div><div> left preconditioning</div><div> using NONE norm type for convergence test</div><div> maximum iterations=2</div><div> tolerances: relative=1e-05, absolute=1e-50, divergence=10000</div><div> left preconditioning</div><div> using nonzero initial guess</div><div> using NONE norm type for convergence test</div><div> PC Object: (mg_levels_4_) 8 MPI processes</div><div> type: sor</div><div> SOR: type = local_symmetric, iterations = 1, local iterations = 1, omega = 1</div><div> linear system matrix = precond matrix:</div><div> Mat Object: () 8 MPI processes</div><div> type: mpiaij</div><div> rows=1532187, cols=1532187, bs=3</div><div> total: nonzeros=1.21304e+08, allocated nonzeros=1.21304e+08</div><div> total number of mallocs used during MatSetValues calls =0</div><div> has attached near null space</div><div> using I-node (on process 0) routines: found 66403 nodes, limit used is 5</div><div> Up solver (post-smoother) same as down solver (pre-smoother)</div><div> linear system matrix = precond matrix:</div><div> Mat Object: () 8 MPI processes</div><div> type: mpiaij</div><div> rows=1532187, cols=1532187, bs=3</div><div> total: nonzeros=1.21304e+08, allocated nonzeros=1.21304e+08</div><div> total number of mallocs used during MatSetValues calls =0</div><div> has attached near null space</div><div> using I-node (on process 0) routines: found 66403 nodes, limit used is 5</div><div>Error, conv_flag < 0!</div></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><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 class="gmail_extra"><div class="gmail_quote"><div><div><div><br></div><div><br></div><div> </div><div><br><div><div><div class="gmail_extra"><br><div class="gmail_quote">On Wed, Nov 11, 2015 at 1:36 PM, David Knezevic <span dir="ltr"><<a href="mailto:david.knezevic@akselos.com" target="_blank">david.knezevic@akselos.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"><div class="gmail_extra"><div class="gmail_quote"><div><div>On Wed, Nov 11, 2015 at 12:57 PM, David Knezevic <span dir="ltr"><<a href="mailto:david.knezevic@akselos.com" target="_blank">david.knezevic@akselos.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"><div><div><div class="gmail_extra"><div class="gmail_quote">On Wed, Nov 11, 2015 at 12:24 PM, David Knezevic <span dir="ltr"><<a href="mailto:david.knezevic@akselos.com" target="_blank">david.knezevic@akselos.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"><div class="gmail_extra"><div><div><div class="gmail_quote">On Tue, Nov 10, 2015 at 10:28 PM, David Knezevic <span dir="ltr"><<a href="mailto:david.knezevic@akselos.com" target="_blank">david.knezevic@akselos.com</a>></span> wrote:<br></div></div></div><div class="gmail_quote"><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 class="gmail_extra"><div class="gmail_quote"><div><div>On Tue, Nov 10, 2015 at 10:24 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: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 class="gmail_extra"><div class="gmail_quote"><div><div>On Tue, Nov 10, 2015 at 9:21 PM, David Knezevic <span dir="ltr"><<a href="mailto:david.knezevic@akselos.com" target="_blank">david.knezevic@akselos.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"><div class="gmail_extra"><div class="gmail_quote">On Tue, Nov 10, 2015 at 10:00 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: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 class="gmail_extra"><div class="gmail_quote"><div><div>On Tue, Nov 10, 2015 at 8:39 PM, David Knezevic <span dir="ltr"><<a href="mailto:david.knezevic@akselos.com" target="_blank">david.knezevic@akselos.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"><div>I'm looking into using GAMG, so I wanted to start with a simple 3D elasticity problem. When I first tried this, I got the following "zero pivot" error:</div><div><br></div><div>-----------------------------------------------------------------------</div><div><br></div><div><div>[0]PETSC ERROR: Zero pivot in LU factorization: <a href="http://www.mcs.anl.gov/petsc/documentation/faq.html#zeropivot" target="_blank">http://www.mcs.anl.gov/petsc/documentation/faq.html#zeropivot</a></div><div>[0]PETSC ERROR: Zero pivot, row 3</div><div>[0]PETSC ERROR: See <a href="http://www.mcs.anl.gov/petsc/documentation/faq.html" target="_blank">http://www.mcs.anl.gov/petsc/documentation/faq.html</a> for trouble shooting.</div><div>[0]PETSC ERROR: Petsc Release Version 3.6.1, Jul, 22, 2015 </div><div>[0]PETSC ERROR: /home/dknez/akselos-dev/scrbe/build/bin/fe_solver-opt_real on a arch-linux2-c-opt named david-Lenovo by dknez Tue Nov 10 21:26:39 2015</div><div>[0]PETSC ERROR: Configure options --with-shared-libraries=1 --with-debugging=0 --download-suitesparse --download-parmetis --download-blacs --with-blas-lapack-dir=/opt/intel/system_studio_2015.2.050/mkl --CXXFLAGS=-Wl,--no-as-needed --download-scalapack --download-mumps --download-metis --download-superlu_dist --prefix=/home/dknez/software/libmesh_install/opt_real/petsc --download-hypre --download-ml</div><div>[0]PETSC ERROR: #1 PetscKernel_A_gets_inverse_A_5() line 48 in /home/dknez/software/petsc-3.6.1/src/mat/impls/baij/seq/dgefa5.c</div><div>[0]PETSC ERROR: #2 MatSOR_SeqAIJ_Inode() line 2808 in /home/dknez/software/petsc-3.6.1/src/mat/impls/aij/seq/inode.c</div><div>[0]PETSC ERROR: #3 MatSOR() line 3697 in /home/dknez/software/petsc-3.6.1/src/mat/interface/matrix.c</div><div>[0]PETSC ERROR: #4 PCApply_SOR() line 37 in /home/dknez/software/petsc-3.6.1/src/ksp/pc/impls/sor/sor.c</div><div>[0]PETSC ERROR: #5 PCApply() line 482 in /home/dknez/software/petsc-3.6.1/src/ksp/pc/interface/precon.c</div><div>[0]PETSC ERROR: #6 KSP_PCApply() line 242 in /home/dknez/software/petsc-3.6.1/include/petsc/private/kspimpl.h</div><div>[0]PETSC ERROR: #7 KSPInitialResidual() line 63 in /home/dknez/software/petsc-3.6.1/src/ksp/ksp/interface/itres.c</div><div>[0]PETSC ERROR: #8 KSPSolve_GMRES() line 235 in /home/dknez/software/petsc-3.6.1/src/ksp/ksp/impls/gmres/gmres.c</div><div>[0]PETSC ERROR: #9 KSPSolve() line 604 in /home/dknez/software/petsc-3.6.1/src/ksp/ksp/interface/itfunc.c</div><div>[0]PETSC ERROR: #10 KSPSolve_Chebyshev() line 381 in /home/dknez/software/petsc-3.6.1/src/ksp/ksp/impls/cheby/cheby.c</div><div>[0]PETSC ERROR: #11 KSPSolve() line 604 in /home/dknez/software/petsc-3.6.1/src/ksp/ksp/interface/itfunc.c</div><div>[0]PETSC ERROR: #12 PCMGMCycle_Private() line 19 in /home/dknez/software/petsc-3.6.1/src/ksp/pc/impls/mg/mg.c</div><div>[0]PETSC ERROR: #13 PCMGMCycle_Private() line 48 in /home/dknez/software/petsc-3.6.1/src/ksp/pc/impls/mg/mg.c</div><div>[0]PETSC ERROR: #14 PCApply_MG() line 338 in /home/dknez/software/petsc-3.6.1/src/ksp/pc/impls/mg/mg.c</div><div>[0]PETSC ERROR: #15 PCApply() line 482 in /home/dknez/software/petsc-3.6.1/src/ksp/pc/interface/precon.c</div><div>[0]PETSC ERROR: #16 KSP_PCApply() line 242 in /home/dknez/software/petsc-3.6.1/include/petsc/private/kspimpl.h</div><div>[0]PETSC ERROR: #17 KSPSolve_CG() line 139 in /home/dknez/software/petsc-3.6.1/src/ksp/ksp/impls/cg/cg.c</div><div>[0]PETSC ERROR: #18 KSPSolve() line 604 in /home/dknez/software/petsc-3.6.1/src/ksp/ksp/interface/itfunc.c</div></div><div><br></div><div>-----------------------------------------------------------------------<br></div><div><br></div><div>I saw that there was a thread about this in September (subject: "gamg and zero pivots"), and that the fix is to use "<span style="color:rgb(0,0,0);white-space:pre-wrap">-mg_levels_pc_type jacobi." When I do that, the solve succeeds (I pasted the -ksp_view at the end of this email).</span></div><div><span style="color:rgb(0,0,0);white-space:pre-wrap"><br></span></div><div><span style="color:rgb(0,0,0);white-space:pre-wrap">So I have two questions about this:</span></div><div><span style="color:rgb(0,0,0);white-space:pre-wrap"><br></span></div><div>1. Is it surprising that I hit this issue for a 3D elasticity problem? Note that matrix assembly was done in libMesh, I can look into the structure of the assembled matrix more carefully, if needed. Also, note that I can solve this problem with direct solvers just fine.</div></div></blockquote><div><br></div></div></div><div>Yes, this seems like a bug, but it could be some strange BC thing I do not understand.</div></div></div></div></blockquote><div><br></div><div><br></div><div>OK, I can look into the matrix in more detail. I agree that it should have a non-zero diagonal, so I'll have a look at what's happening with that.</div><div><br></div><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 class="gmail_extra"><div class="gmail_quote"><div>Naively, the elastic element matrix has a nonzero diagonal. I see that you are doing LU<br></div><div>of size 5. That seems strange for 3D elasticity. Am I missing something? I would expect</div><div>block size 3.</div></div></div></div></blockquote><div><br></div><div><br></div><div>I'm not sure what is causing the LU of size 5. Is there a setting to control that?<br></div><div><br></div><div>Regarding the block size: I set the vector and matrix block size to 3 via VecSetBlockSize and MatSetBlockSize. I also used MatNullSpaceCreateRigidBody on a vector with block size of 3, and set the matrix's near nullspace using that.</div></div></div></div></blockquote><div><br></div></div></div><div>Can you run this same example with -mat_no_inode? I think it may be a strange blocking that is causing this.</div></div></div></div></blockquote><div><br></div><div><br></div></div></div><div>That works. The -ksp_view output is below.</div></div></div></div></blockquote><div><br></div><div><br></div></div></div><div>I just wanted to follow up on this. I had a more careful look at the matrix, and confirmed that there are no zero entries on the diagonal (as expected for elasticity). The matrix is from one of libMesh's example problems: a simple cantilever model using HEX8 elements.</div><div><br></div><div>Do you have any further thoughts about what might cause the "strange blocking" that you referred to? If there's something non-standard that libMesh is doing with the blocks, I'd be interested to look into that. I can send over the matrix if that would be helpful.</div><div><br></div><div>Thanks,<br>David</div><div><br></div></div></div></div>
</blockquote></div><br></div></div></div><div class="gmail_extra">P.S. I was previously calling VecSetBlockSize and MatSetBlockSize to set the block size to 3. When I don't do that, I no longer need to call -mat_no_inodes. I've pasted the -ksp_view output below. Does it look like that's working OK?</div></div></blockquote><div><br></div><div><br></div></div></div><div>Sorry for the multiple messages, but I think I found the issue. libMesh internally sets the block size to 1 earlier on (in PetscMatrix::init()). I guess it'll work fine if I get it to set the block size to 3 instead, so I'll look into that. (libMesh has an enable-blocked-storage configure option that should take care of this automatically.)<span><font color="#888888"><br></font></span></div><span><font color="#888888"><div><br></div><div>David</div><div><br></div><div><br></div></font></span></div></div></div>
</blockquote></div><br></div>
</div></div></div></div></div></div><br></div></div>
</blockquote></div><br></div></div>
</blockquote></div><br><br clear="all"><div><br></div>-- <br><div class="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>