On Fri, Nov 13, 2015 at 9:45 AM, Matthew Knepley <[email protected]> wrote:
> On Thu, Nov 12, 2015 at 10:02 PM, David Knezevic < > [email protected]> wrote: > >> On Thu, Nov 12, 2015 at 9:58 AM, David Knezevic < >> [email protected]> wrote: >> >>> On Thu, Nov 12, 2015 at 9:50 AM, Mark Adams <[email protected]> wrote: >>> >>>> 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. >>>> >>>> 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. >>>> >>> >>> OK, thanks for the info, I'll look into this further. >>> >>> 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. >>> >> >> >> As I mentioned above, GAMG seems to be working fine for me now. Thanks >> for the help on this. >> >> 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.) >> >> Thanks, >> David >> >> -------------------------------------------- >> >> 0 KSP Residual norm 2.475892877233e-03 >> 1 KSP Residual norm 6.181785698923e-05 >> 2 KSP Residual norm 9.439050821656e-04 >> > > Something very strange is happening here. CG should converge > monotonically, but above it does not. What could be happening? > > a) We could have lost orthogonality > > This seems really unlikely in 2 iterates. > > b) The preconditioner could be nonlinear > > In order to test this, could you run with both GMRES and FGMRES? If > the convergence is different, then something is wrong > with the preconditioner. > > I have looked below, and it seems like this should be linear. Mark, > could something be happening in GAMG? > > 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 > 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 > get a sense of discretization error using a manufactured solution. > OK, thanks for this info. I'll do some tests based on your comments above when I get some time and I'll let you know what happens. Thanks, David > KSP Object: 8 MPI processes >> type: cg >> maximum iterations=5000 >> tolerances: relative=1e-12, absolute=1e-50, divergence=10000 >> left preconditioning >> using nonzero initial guess >> using PRECONDITIONED norm type for convergence test >> PC Object: 8 MPI processes >> type: gamg >> MG: type is MULTIPLICATIVE, levels=5 cycles=v >> Cycles per PCApply=1 >> Using Galerkin computed coarse grid matrices >> GAMG specific options >> Threshold for dropping small values from graph 0 >> AGG specific options >> Symmetric graph false >> Coarse grid solver -- level ------------------------------- >> KSP Object: (mg_coarse_) 8 MPI processes >> type: gmres >> GMRES: restart=30, using Classical (unmodified) Gram-Schmidt >> Orthogonalization with no iterative refinement >> GMRES: happy breakdown tolerance 1e-30 >> maximum iterations=1, initial guess is zero >> tolerances: relative=1e-05, absolute=1e-50, divergence=10000 >> left preconditioning >> using NONE norm type for convergence test >> PC Object: (mg_coarse_) 8 MPI processes >> type: bjacobi >> block Jacobi: number of blocks = 8 >> Local solve is same for all blocks, in the following KSP and PC >> objects: >> KSP Object: (mg_coarse_sub_) 1 MPI processes >> type: preonly >> maximum iterations=1, initial guess is zero >> tolerances: relative=1e-05, absolute=1e-50, divergence=10000 >> left preconditioning >> using NONE norm type for convergence test >> PC Object: (mg_coarse_sub_) 1 MPI processes >> type: lu >> LU: out-of-place factorization >> tolerance for zero pivot 2.22045e-14 >> using diagonal shift on blocks to prevent zero pivot [INBLOCKS] >> matrix ordering: nd >> factor fill ratio given 5, needed 1 >> Factored matrix follows: >> Mat Object: 1 MPI processes >> type: seqaij >> rows=6, cols=6, bs=6 >> package used to perform factorization: petsc >> total: nonzeros=36, allocated nonzeros=36 >> total number of mallocs used during MatSetValues calls =0 >> using I-node routines: found 2 nodes, limit used is 5 >> linear system matrix = precond matrix: >> Mat Object: 1 MPI processes >> type: seqaij >> rows=6, cols=6, bs=6 >> total: nonzeros=36, allocated nonzeros=36 >> total number of mallocs used during MatSetValues calls =0 >> using I-node routines: found 2 nodes, limit used is 5 >> linear system matrix = precond matrix: >> Mat Object: 8 MPI processes >> type: mpiaij >> rows=6, cols=6, bs=6 >> total: nonzeros=36, allocated nonzeros=36 >> total number of mallocs used during MatSetValues calls =0 >> using I-node (on process 0) routines: found 2 nodes, limit used >> is 5 >> Down solver (pre-smoother) on level 1 ------------------------------- >> KSP Object: (mg_levels_1_) 8 MPI processes >> type: chebyshev >> Chebyshev: eigenvalue estimates: min = 0.146922, max = 1.61615 >> Chebyshev: eigenvalues estimated using gmres with translations >> [0 0.1; 0 1.1] >> KSP Object: (mg_levels_1_esteig_) 8 MPI processes >> type: gmres >> GMRES: restart=30, using Classical (unmodified) Gram-Schmidt >> Orthogonalization with no iterative refinement >> GMRES: happy breakdown tolerance 1e-30 >> maximum iterations=10, initial guess is zero >> tolerances: relative=1e-05, absolute=1e-50, divergence=10000 >> left preconditioning >> using NONE norm type for convergence test >> maximum iterations=2 >> tolerances: relative=1e-05, absolute=1e-50, divergence=10000 >> left preconditioning >> using nonzero initial guess >> using NONE norm type for convergence test >> PC Object: (mg_levels_1_) 8 MPI processes >> type: sor >> SOR: type = local_symmetric, iterations = 1, local iterations = >> 1, omega = 1 >> linear system matrix = precond matrix: >> Mat Object: 8 MPI processes >> type: mpiaij >> rows=162, cols=162, bs=6 >> total: nonzeros=26244, allocated nonzeros=26244 >> total number of mallocs used during MatSetValues calls =0 >> using I-node (on process 0) routines: found 5 nodes, limit used >> is 5 >> Up solver (post-smoother) same as down solver (pre-smoother) >> Down solver (pre-smoother) on level 2 ------------------------------- >> KSP Object: (mg_levels_2_) 8 MPI processes >> type: chebyshev >> Chebyshev: eigenvalue estimates: min = 0.143941, max = 1.58335 >> Chebyshev: eigenvalues estimated using gmres with translations >> [0 0.1; 0 1.1] >> KSP Object: (mg_levels_2_esteig_) 8 MPI processes >> type: gmres >> GMRES: restart=30, using Classical (unmodified) Gram-Schmidt >> Orthogonalization with no iterative refinement >> GMRES: happy breakdown tolerance 1e-30 >> maximum iterations=10, initial guess is zero >> tolerances: relative=1e-05, absolute=1e-50, divergence=10000 >> left preconditioning >> using NONE norm type for convergence test >> maximum iterations=2 >> tolerances: relative=1e-05, absolute=1e-50, divergence=10000 >> left preconditioning >> using nonzero initial guess >> using NONE norm type for convergence test >> PC Object: (mg_levels_2_) 8 MPI processes >> type: sor >> SOR: type = local_symmetric, iterations = 1, local iterations = >> 1, omega = 1 >> linear system matrix = precond matrix: >> Mat Object: 8 MPI processes >> type: mpiaij >> rows=4668, cols=4668, bs=6 >> total: nonzeros=3.00254e+06, allocated nonzeros=3.00254e+06 >> total number of mallocs used during MatSetValues calls =0 >> using I-node (on process 0) routines: found 158 nodes, limit >> used is 5 >> Up solver (post-smoother) same as down solver (pre-smoother) >> Down solver (pre-smoother) on level 3 ------------------------------- >> KSP Object: (mg_levels_3_) 8 MPI processes >> type: chebyshev >> Chebyshev: eigenvalue estimates: min = 0.128239, max = 1.41063 >> Chebyshev: eigenvalues estimated using gmres with translations >> [0 0.1; 0 1.1] >> KSP Object: (mg_levels_3_esteig_) 8 MPI processes >> type: gmres >> GMRES: restart=30, using Classical (unmodified) Gram-Schmidt >> Orthogonalization with no iterative refinement >> GMRES: happy breakdown tolerance 1e-30 >> maximum iterations=10, initial guess is zero >> tolerances: relative=1e-05, absolute=1e-50, divergence=10000 >> left preconditioning >> using NONE norm type for convergence test >> maximum iterations=2 >> tolerances: relative=1e-05, absolute=1e-50, divergence=10000 >> left preconditioning >> using nonzero initial guess >> using NONE norm type for convergence test >> PC Object: (mg_levels_3_) 8 MPI processes >> type: sor >> SOR: type = local_symmetric, iterations = 1, local iterations = >> 1, omega = 1 >> linear system matrix = precond matrix: >> Mat Object: 8 MPI processes >> type: mpiaij >> rows=66390, cols=66390, bs=6 >> total: nonzeros=1.65982e+07, allocated nonzeros=1.65982e+07 >> total number of mallocs used during MatSetValues calls =0 >> using I-node (on process 0) routines: found 2668 nodes, limit >> used is 5 >> Up solver (post-smoother) same as down solver (pre-smoother) >> Down solver (pre-smoother) on level 4 ------------------------------- >> KSP Object: (mg_levels_4_) 8 MPI processes >> type: chebyshev >> Chebyshev: eigenvalue estimates: min = 0.1, max = 1.1 >> Chebyshev: eigenvalues estimated using gmres with translations >> [0 0.1; 0 1.1] >> KSP Object: (mg_levels_4_esteig_) 8 MPI processes >> type: gmres >> GMRES: restart=30, using Classical (unmodified) Gram-Schmidt >> Orthogonalization with no iterative refinement >> GMRES: happy breakdown tolerance 1e-30 >> maximum iterations=10, initial guess is zero >> tolerances: relative=1e-05, absolute=1e-50, divergence=10000 >> left preconditioning >> using NONE norm type for convergence test >> maximum iterations=2 >> tolerances: relative=1e-05, absolute=1e-50, divergence=10000 >> left preconditioning >> using nonzero initial guess >> using NONE norm type for convergence test >> PC Object: (mg_levels_4_) 8 MPI processes >> type: sor >> SOR: type = local_symmetric, iterations = 1, local iterations = >> 1, omega = 1 >> linear system matrix = precond matrix: >> Mat Object: () 8 MPI processes >> type: mpiaij >> rows=1532187, cols=1532187, bs=3 >> total: nonzeros=1.21304e+08, allocated nonzeros=1.21304e+08 >> total number of mallocs used during MatSetValues calls =0 >> has attached near null space >> using I-node (on process 0) routines: found 66403 nodes, limit >> used is 5 >> Up solver (post-smoother) same as down solver (pre-smoother) >> linear system matrix = precond matrix: >> Mat Object: () 8 MPI processes >> type: mpiaij >> rows=1532187, cols=1532187, bs=3 >> total: nonzeros=1.21304e+08, allocated nonzeros=1.21304e+08 >> total number of mallocs used during MatSetValues calls =0 >> has attached near null space >> using I-node (on process 0) routines: found 66403 nodes, limit used >> is 5 >> Error, conv_flag < 0! >> >> >> >> >> >> >> >> >>> >>> >>> >>> >>> >>> On Wed, Nov 11, 2015 at 1:36 PM, David Knezevic < >>> [email protected]> wrote: >>> >>>> On Wed, Nov 11, 2015 at 12:57 PM, David Knezevic < >>>> [email protected]> wrote: >>>> >>>>> On Wed, Nov 11, 2015 at 12:24 PM, David Knezevic < >>>>> [email protected]> wrote: >>>>> >>>>>> On Tue, Nov 10, 2015 at 10:28 PM, David Knezevic < >>>>>> [email protected]> wrote: >>>>>> >>>>>>> On Tue, Nov 10, 2015 at 10:24 PM, Matthew Knepley <[email protected] >>>>>>> > wrote: >>>>>>> >>>>>>>> On Tue, Nov 10, 2015 at 9:21 PM, David Knezevic < >>>>>>>> [email protected]> wrote: >>>>>>>> >>>>>>>>> On Tue, Nov 10, 2015 at 10:00 PM, Matthew Knepley < >>>>>>>>> [email protected]> wrote: >>>>>>>>> >>>>>>>>>> On Tue, Nov 10, 2015 at 8:39 PM, David Knezevic < >>>>>>>>>> [email protected]> wrote: >>>>>>>>>> >>>>>>>>>>> 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: >>>>>>>>>>> >>>>>>>>>>> >>>>>>>>>>> ----------------------------------------------------------------------- >>>>>>>>>>> >>>>>>>>>>> [0]PETSC ERROR: Zero pivot in LU factorization: >>>>>>>>>>> http://www.mcs.anl.gov/petsc/documentation/faq.html#zeropivot >>>>>>>>>>> [0]PETSC ERROR: Zero pivot, row 3 >>>>>>>>>>> [0]PETSC ERROR: See >>>>>>>>>>> http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble >>>>>>>>>>> shooting. >>>>>>>>>>> [0]PETSC ERROR: Petsc Release Version 3.6.1, Jul, 22, 2015 >>>>>>>>>>> [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 >>>>>>>>>>> [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 >>>>>>>>>>> [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 >>>>>>>>>>> [0]PETSC ERROR: #2 MatSOR_SeqAIJ_Inode() line 2808 in >>>>>>>>>>> /home/dknez/software/petsc-3.6.1/src/mat/impls/aij/seq/inode.c >>>>>>>>>>> [0]PETSC ERROR: #3 MatSOR() line 3697 in >>>>>>>>>>> /home/dknez/software/petsc-3.6.1/src/mat/interface/matrix.c >>>>>>>>>>> [0]PETSC ERROR: #4 PCApply_SOR() line 37 in >>>>>>>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/pc/impls/sor/sor.c >>>>>>>>>>> [0]PETSC ERROR: #5 PCApply() line 482 in >>>>>>>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/pc/interface/precon.c >>>>>>>>>>> [0]PETSC ERROR: #6 KSP_PCApply() line 242 in >>>>>>>>>>> /home/dknez/software/petsc-3.6.1/include/petsc/private/kspimpl.h >>>>>>>>>>> [0]PETSC ERROR: #7 KSPInitialResidual() line 63 in >>>>>>>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/ksp/interface/itres.c >>>>>>>>>>> [0]PETSC ERROR: #8 KSPSolve_GMRES() line 235 in >>>>>>>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/ksp/impls/gmres/gmres.c >>>>>>>>>>> [0]PETSC ERROR: #9 KSPSolve() line 604 in >>>>>>>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/ksp/interface/itfunc.c >>>>>>>>>>> [0]PETSC ERROR: #10 KSPSolve_Chebyshev() line 381 in >>>>>>>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/ksp/impls/cheby/cheby.c >>>>>>>>>>> [0]PETSC ERROR: #11 KSPSolve() line 604 in >>>>>>>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/ksp/interface/itfunc.c >>>>>>>>>>> [0]PETSC ERROR: #12 PCMGMCycle_Private() line 19 in >>>>>>>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/pc/impls/mg/mg.c >>>>>>>>>>> [0]PETSC ERROR: #13 PCMGMCycle_Private() line 48 in >>>>>>>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/pc/impls/mg/mg.c >>>>>>>>>>> [0]PETSC ERROR: #14 PCApply_MG() line 338 in >>>>>>>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/pc/impls/mg/mg.c >>>>>>>>>>> [0]PETSC ERROR: #15 PCApply() line 482 in >>>>>>>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/pc/interface/precon.c >>>>>>>>>>> [0]PETSC ERROR: #16 KSP_PCApply() line 242 in >>>>>>>>>>> /home/dknez/software/petsc-3.6.1/include/petsc/private/kspimpl.h >>>>>>>>>>> [0]PETSC ERROR: #17 KSPSolve_CG() line 139 in >>>>>>>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/ksp/impls/cg/cg.c >>>>>>>>>>> [0]PETSC ERROR: #18 KSPSolve() line 604 in >>>>>>>>>>> /home/dknez/software/petsc-3.6.1/src/ksp/ksp/interface/itfunc.c >>>>>>>>>>> >>>>>>>>>>> >>>>>>>>>>> ----------------------------------------------------------------------- >>>>>>>>>>> >>>>>>>>>>> I saw that there was a thread about this in September (subject: >>>>>>>>>>> "gamg and zero pivots"), and that the fix is to use >>>>>>>>>>> "-mg_levels_pc_type >>>>>>>>>>> jacobi." When I do that, the solve succeeds (I pasted the -ksp_view >>>>>>>>>>> at the >>>>>>>>>>> end of this email). >>>>>>>>>>> >>>>>>>>>>> So I have two questions about this: >>>>>>>>>>> >>>>>>>>>>> 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. >>>>>>>>>>> >>>>>>>>>> >>>>>>>>>> Yes, this seems like a bug, but it could be some strange BC thing >>>>>>>>>> I do not understand. >>>>>>>>>> >>>>>>>>> >>>>>>>>> >>>>>>>>> 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. >>>>>>>>> >>>>>>>>> >>>>>>>>> >>>>>>>>> >>>>>>>>>> Naively, the elastic element matrix has a nonzero diagonal. I see >>>>>>>>>> that you are doing LU >>>>>>>>>> of size 5. That seems strange for 3D elasticity. Am I missing >>>>>>>>>> something? I would expect >>>>>>>>>> block size 3. >>>>>>>>>> >>>>>>>>> >>>>>>>>> >>>>>>>>> I'm not sure what is causing the LU of size 5. Is there a setting >>>>>>>>> to control that? >>>>>>>>> >>>>>>>>> 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. >>>>>>>>> >>>>>>>> >>>>>>>> Can you run this same example with -mat_no_inode? I think it may be >>>>>>>> a strange blocking that is causing this. >>>>>>>> >>>>>>> >>>>>>> >>>>>>> That works. The -ksp_view output is below. >>>>>>> >>>>>> >>>>>> >>>>>> 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. >>>>>> >>>>>> 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. >>>>>> >>>>>> Thanks, >>>>>> David >>>>>> >>>>>> >>>>> 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? >>>>> >>>> >>>> >>>> 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.) >>>> >>>> David >>>> >>>> >>>> >>> >>> >> > > > -- > What most experimenters take for granted before they begin their > experiments is infinitely more interesting than any results to which their > experiments lead. > -- Norbert Wiener >
