<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">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>Can you run this same example with -mat_no_inode? I think it may be a strange blocking that is causing this.</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-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"><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><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>2. Is there a way to set "<span style="color:rgb(0,0,0);white-space:pre-wrap">-mg_levels_pc_type jacobi" programmatically, rather than via the command line?</span></div></div></blockquote><div><br></div></span><div>I would really discourage you from doing this. It makes your code fragile and inflexible.</div></div></div></div></blockquote><div><br></div><div>OK. The reason I asked is that in this case I have to write a bunch of code to set block sizes and the near nullspace, so I figured it'd be good to also set the corresponding solver options required to make this work... anyway, if I can fix the zero diagonal issue, I guess this will be moot.</div><div><br></div><div>David</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"><span><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>-----------------------------------------------------------------------<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"><br></span></div><div><span style="color:rgb(0,0,0);white-space:pre-wrap">ksp_view output:</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"><br></span></div><div><font color="#000000"><span style="white-space:pre-wrap">KSP Object: 1 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: 1 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_)     1 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_)     1 MPI processes
      type: bjacobi
        block Jacobi: number of blocks = 1
        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=30, cols=30, bs=6
                  package used to perform factorization: petsc
                  total: nonzeros=540, allocated nonzeros=540
                  total number of mallocs used during MatSetValues calls =0
                    using I-node routines: found 9 nodes, limit used is 5
          linear system matrix = precond matrix:
          Mat Object:           1 MPI processes
            type: seqaij
            rows=30, cols=30, bs=6
            total: nonzeros=540, allocated nonzeros=540
            total number of mallocs used during MatSetValues calls =0
              using I-node routines: found 9 nodes, limit used is 5
      linear system matrix = precond matrix:
      Mat Object:       1 MPI processes
        type: seqaij
        rows=30, cols=30, bs=6
        total: nonzeros=540, allocated nonzeros=540
        total number of mallocs used during MatSetValues calls =0
          using I-node routines: found 9 nodes, limit used is 5
  Down solver (pre-smoother) on level 1 -------------------------------
    KSP Object:    (mg_levels_1_)     1 MPI processes
      type: chebyshev
        Chebyshev: eigenvalue estimates:  min = 0.335276, max = 3.68804
        Chebyshev: eigenvalues estimated using gmres with translations  [0 0.1; 0 1.1]
        KSP Object:        (mg_levels_1_esteig_)         1 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_)     1 MPI processes
      type: jacobi
      linear system matrix = precond matrix:
      Mat Object:       1 MPI processes
        type: seqaij
        rows=72, cols=72, bs=6
        total: nonzeros=1728, allocated nonzeros=1728
        total number of mallocs used during MatSetValues calls =0
          using I-node routines: found 23 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_)     1 MPI processes
      type: chebyshev
        Chebyshev: eigenvalue estimates:  min = 0.260121, max = 2.86133
        Chebyshev: eigenvalues estimated using gmres with translations  [0 0.1; 0 1.1]
        KSP Object:        (mg_levels_2_esteig_)         1 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_)     1 MPI processes
      type: jacobi
      linear system matrix = precond matrix:
      Mat Object:       1 MPI processes
        type: seqaij
        rows=174, cols=174, bs=6
        total: nonzeros=5796, allocated nonzeros=5796
        total number of mallocs used during MatSetValues calls =0
          using I-node routines: found 57 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_)     1 MPI processes
      type: chebyshev
        Chebyshev: eigenvalue estimates:  min = 0.267401, max = 2.94141
        Chebyshev: eigenvalues estimated using gmres with translations  [0 0.1; 0 1.1]
        KSP Object:        (mg_levels_3_esteig_)         1 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_)     1 MPI processes
      type: jacobi
      linear system matrix = precond matrix:
      Mat Object:       1 MPI processes
        type: seqaij
        rows=828, cols=828, bs=6
        total: nonzeros=44496, allocated nonzeros=44496
        total number of mallocs used during MatSetValues calls =0
          using I-node routines: found 276 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_)     1 MPI processes
      type: chebyshev
        Chebyshev: eigenvalue estimates:  min = 0.224361, max = 2.46797
        Chebyshev: eigenvalues estimated using gmres with translations  [0 0.1; 0 1.1]
        KSP Object:        (mg_levels_4_esteig_)         1 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_)     1 MPI processes
      type: jacobi
      linear system matrix = precond matrix:
      Mat Object:      ()       1 MPI processes
        type: seqaij
        rows=2676, cols=2676, bs=3
        total: nonzeros=94014, allocated nonzeros=94014
        total number of mallocs used during MatSetValues calls =0
          has attached near null space
          using I-node routines: found 892 nodes, limit used is 5
  Up solver (post-smoother) same as down solver (pre-smoother)
  linear system matrix = precond matrix:
  Mat Object:  ()   1 MPI processes
    type: seqaij
    rows=2676, cols=2676, bs=3
    total: nonzeros=94014, allocated nonzeros=94014
    total number of mallocs used during MatSetValues calls =0
      has attached near null space
      using I-node routines: found 892 nodes, limit used is 5</span></font><span style="color:rgb(0,0,0);white-space:pre-wrap">
</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"><br></span></div></div>
</blockquote></span></div><span><font color="#888888"><br><br clear="all"><span class=""><font color="#888888"><div><br></div>-- <br><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>
</font></span></font></span></div></div>
</blockquote></div></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>