get some info of the solver status of PETSc

MathLib/LinAlg/PETSc/PETScLinearSolver.cpp

@@ -25,15 +25,15 @@ PETScLinearSolver::PETScLinearSolver(PETScMatrix &A, const std::string &prefix)
     KSPGetPC(_solver, &_pc);
 
     //
-    if ( !prefix.empty() ) 
-    {       
-       KSPSetOptionsPrefix(_solver, prefix.c_str());
+    if ( !prefix.empty() )
+    {
+        KSPSetOptionsPrefix(_solver, prefix.c_str());
     }
-               
+
     KSPSetFromOptions(_solver);  // set running time option
 }
 
-void PETScLinearSolver::solve(const PETScVector &b, PETScVector &x)
+bool PETScLinearSolver::solve(const PETScVector &b, PETScVector &x)
 {
 // define TEST_MEM_PETSC
 #ifdef TEST_MEM_PETSC
@@ -48,16 +48,8 @@ void PETScLinearSolver::solve(const PETScVector &b, PETScVector &x)
     KSPConvergedReason reason;
     KSPGetConvergedReason(_solver, &reason);
 
-    if(reason == KSP_DIVERGED_INDEFINITE_PC)
-    {
-        PetscPrintf(PETSC_COMM_WORLD,"\nDivergence because of indefinite preconditioner;\n");
-        PetscPrintf(PETSC_COMM_WORLD,"Run the executable again but with -pc_factor_shift_positive_definite option.\n");
-    }
-    else if(reason < 0)
-    {
-        PetscPrintf(PETSC_COMM_WORLD,"\nOther kind of divergence: this should not happen.\n");
-    }
-    else
+    bool converged = true;
+    if(reason > 0)
     {
         const char *slv_type;
         const char *prc_type;
@@ -71,14 +63,47 @@ void PETScLinearSolver::solve(const PETScVector &b, PETScVector &x)
         int its;
         KSPGetIterationNumber(_solver, &its);
         PetscPrintf(PETSC_COMM_WORLD,"\nConvergence in %d iterations.\n", its);
-        PetscPrintf(PETSC_COMM_WORLD,"\n================================================");
+        PetscPrintf(PETSC_COMM_WORLD,"\n================================================\n");
+    }
+    else if(reason == KSP_DIVERGED_ITS)
+    {
+        PetscPrintf(PETSC_COMM_WORLD, "\nWaning: maximum iteration reached.\n");
+        converged = false;
+    }
+    else
+    {
+        if(reason == KSP_DIVERGED_INDEFINITE_PC)
+        {
+            PetscPrintf(PETSC_COMM_WORLD, "\nDivergence because of indefinite preconditioner,");
+            PetscPrintf(PETSC_COMM_WORLD, "\nTry to run again with -pc_factor_shift_positive_definite option.\n");
+        }
+        else if(reason == KSP_DIVERGED_BREAKDOWN_BICG)
+        {
+            PetscPrintf(PETSC_COMM_WORLD, "\nKSPBICG method was detected so the method could not continue to enlarge the Krylov space.");
+            PetscPrintf(PETSC_COMM_WORLD, "\nTry to run again with another solver.\n");
+        }
+
+        else if(reason == KSP_DIVERGED_NONSYMMETRIC)
+        {
+            PetscPrintf(PETSC_COMM_WORLD, "\nMatrx or preconditioner is unsymmetric but KSP requires symmetric.\n");
+        }
+        else
+        {
+            PetscPrintf(PETSC_COMM_WORLD, "\nOther kind divergence, use command option -ksp_monitor or -log_summary to check the details.\n");
+        }
+
+        PetscPrintf(PETSC_COMM_WORLD, "\nLinear solver (PETSc KSP) failed, quit now.\n");
+        KSPDestroy(&_solver);
+        PetscFinalize();
+        exit(EXIT_FAILURE);
     }
-    PetscPrintf(PETSC_COMM_WORLD,"\n");
 
 #ifdef TEST_MEM_PETSC
     PetscMemoryGetCurrentUsage(&mem2);
     PetscPrintf(PETSC_COMM_WORLD, "###Memory usage by solver. Before :%f After:%f Increase:%d\n", mem1, mem2, (int)(mem2 - mem1));
 #endif
+
+    return converged;
 }
 
 } //end of namespace

MathLib/LinAlg/PETSc/PETScLinearSolver.h

@@ -54,16 +54,30 @@ class PETScLinearSolver
             Solve a system of equations.
             \param b The right hand of the equations.
             \param x The solutions to be solved.
+            \return  true: converged, false: diverged due to exeeding
+                     the maximum iterations.
         */
-        void solve(const PETScVector &b, PETScVector &x);
+        bool solve(const PETScVector &b, PETScVector &x);
+
+        /*!
+            \brief Get number of iterations.
+            If function solve(...) returns false, it returns the maximum
+            iterations.
+        */
+        PetscInt getIterations() const
+        {
+            PetscInt its = 0;
+            KSPGetIterationNumber(_solver, &its);
+            return its;
+        }
 
     private:
         /// Matrix, kept as a member only for solving successive linear equation
-        /// that preconditioner matrix may vary.    
-        PETScMatrix &_A;  
-    
+        /// that preconditioner matrix may vary.
+        PETScMatrix &_A;
+
         KSP _solver; ///< Slover type.
-        PC _pc;      ///< Preconditioner type.        
+        PC _pc;      ///< Preconditioner type.
 };
 
 } // end namespace

MathLib/LinAlg/PETSc/PETScTools.cpp

@@ -28,10 +28,10 @@ void applyKnownSolution(PETScMatrix &A, PETScVector &b, PETScVector &x,
 
     if(vec_knownX_id.size() > 0)
     {
-       x.set(vec_knownX_id, vec_knownX_x);
-       b.set(vec_knownX_id, vec_knownX_x);
+        x.set(vec_knownX_id, vec_knownX_x);
+        b.set(vec_knownX_id, vec_knownX_x);
     }
-     
+
     x.finalizeAssembly();
     b.finalizeAssembly();
 }

MathLib/LinAlg/PETSc/PETScTools.h

@@ -26,7 +26,7 @@ namespace MathLib
 {
 /*!
    \brief apply known solutions to a system of linear equations
-  
+
    \param A                 Coefficient matrix
    \param b                 RHS vector
    \param vec_knownX_id    a vector of known solution entry IDs

Tests/MathLib/TestLinearSolver.cpp

@@ -155,7 +155,7 @@ void checkLinearSolverInterface(T_MATRIX &A, T_VECTOR &b, const std::string &pre
     std::vector<double> local_vec(2);
     local_vec[0] = mrank+1;
     local_vec[1] = 2. * (mrank+1);
-    x.set(col_pos, local_vec);
+    x.set(row_pos, local_vec);
 
     double x0[6];
     double x1[6];
@@ -164,7 +164,7 @@ void checkLinearSolverInterface(T_MATRIX &A, T_VECTOR &b, const std::string &pre
     A.multiply(x, b);
 
     // apply BC
-    std::vector<int> bc_id;  /// Type must be int to match Petsc_Int
+    std::vector<int> bc_id;  // Type must be int to match Petsc_Int
     std::vector<double> bc_value;
 
     if(mrank == 1)
@@ -182,9 +182,11 @@ void checkLinearSolverInterface(T_MATRIX &A, T_VECTOR &b, const std::string &pre
     // solve
     T_LINEAR_SOVLER ls(A, prefix_name);
     ls.solve(b, x);
+    
+    EXPECT_GT(ls.getIterations(), 0u);
+    
     x.getGlobalVector(x1);
-
-    ASSERT_ARRAY_NEAR(x0, x1, 6, 1e-5);
+    ASSERT_ARRAY_NEAR(x0, x1, 6, 1e-5);        
 }
 #endif