Added options to visualise matrix

diff --git a/CoreFlows/Models/inc/ProblemCoreFlows.hxx b/CoreFlows/Models/inc/ProblemCoreFlows.hxx
index ee2245e6308230a4b484444e59f5c5a6193ac635..738f0496360bbcd0a0de776fb765baa270ea245e 100755 (executable)
--- a/CoreFlows/Models/inc/ProblemCoreFlows.hxx
+++ b/CoreFlows/Models/inc/ProblemCoreFlows.hxx
@@ -731,10 +731,10 @@ public :
 
     //Spectral analysis
     double getConditionNumber(bool isSingular=false, double tol=1e-6) const;
-    std::vector< double > getEigenvalues (int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const;
+    std::vector< double > getEigenvalues (int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6, EPSType type = EPSKRYLOVSCHUR, bool viewEigenvaluesInXWindows=false, double pause_lenght=0) const;
     std::vector< Vector > getEigenvectors(int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const;
     Field getEigenvectorsField(int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const;
-    std::vector< double > getSingularValues( int nsv, SVDWhich which=SVD_SMALLEST, double tol=1e-6) const;
+    std::vector< double > getSingularValues( int nsv, SVDWhich which=SVD_SMALLEST, double tol=1e-6, SVDType type = SVDCYCLIC, bool viewSingularValuesInXWindows=false, double pause_lenght=0) const;
     std::vector< Vector > getSingularVectors(int nsv, SVDWhich which=SVD_SMALLEST, double tol=1e-6) const;
 
        //  some supplementary functions
@@ -746,7 +746,7 @@ public :
         * @param name, string, name or description of the matrix
         * @return displays the matrix on the terminal
         *  */
-       static void displayMatrix(double *matrix, int size, string name="Vector coefficients :");
+       static void displayMatrix(double *matrix, int size, string name="Matrix coefficients :");
 
        /** \fn displayMatrix
         * \brief displays a vector of size "size" for profiling
@@ -755,10 +755,11 @@ public :
         * @param name, string, name or description of the vector
         * @return displays the vector on the terminal
         *  */
-       static void displayVector(double *vector, int size, string name="Matrix coefficients :");
+       static void displayVector(double *vector, int size, string name="Vector coefficients :");
 
 protected :
 
+       // Mesh info
        int _Ndim;//space dimension
        int _nVar;//Number of equations to sole
        int _Nmailles;//number of cells
@@ -802,7 +803,7 @@ protected :
         *  */
        void createKSP();
 
-       //simulation monitoring variables
+       // Simulation monitoring variables
        bool _isStationary;
        bool _initialDataSet;
        bool _initializedMemory;

diff --git a/CoreFlows/Models/src/ProblemCoreFlows.cxx b/CoreFlows/Models/src/ProblemCoreFlows.cxx
index 6c35c630964c6f83e8b4f4d13fb21f802844a895..492dc34c648025d9de07f72798e66e212389b3d2 100755 (executable)
--- a/CoreFlows/Models/src/ProblemCoreFlows.cxx
+++ b/CoreFlows/Models/src/ProblemCoreFlows.cxx
@@ -732,10 +732,10 @@ ProblemCoreFlows::getConditionNumber(bool isSingular, double tol) const
   return A.getConditionNumber( isSingular, tol);
 }
 std::vector< double > 
-ProblemCoreFlows::getEigenvalues(int nev, EPSWhich which, double tol) const
+ProblemCoreFlows::getEigenvalues(int nev, EPSWhich which, double tol, EPSType type, bool viewEigenvaluesInXWindows, double pause_lenght) const
 {
   SparseMatrixPetsc A = SparseMatrixPetsc(_A);
-  return A.getEigenvalues( nev, which, tol);
+  return A.getEigenvalues( nev, which, tol, type, viewEigenvaluesInXWindows, pause_lenght);
 }
 std::vector< Vector > 
 ProblemCoreFlows::getEigenvectors(int nev, EPSWhich which, double tol) const
@@ -758,10 +758,10 @@ ProblemCoreFlows::getEigenvectorsField(int nev, EPSWhich which, double tol) cons
 }
 
 std::vector< double > 
-ProblemCoreFlows::getSingularValues( int nsv, SVDWhich which, double tol) const
+ProblemCoreFlows::getSingularValues( int nsv, SVDWhich which, double tol, SVDType type, bool viewSingularValuesInXWindows, double pause_lenght) const
 {
   SparseMatrixPetsc A = SparseMatrixPetsc(_A);
-  return A.getSingularValues( nsv, which, tol);
+  return A.getSingularValues( nsv, which, tol, type, viewSingularValuesInXWindows, pause_lenght);
 }
 std::vector< Vector > 
 ProblemCoreFlows::getSingularVectors(int nsv, SVDWhich which, double tol) const
@@ -810,3 +810,42 @@ ProblemCoreFlows::setHeatPowerField(string fileName, string fieldName, int itera
        _heatPowerField.getMesh().checkFastEquivalWith(_mesh);
        _heatPowerFieldSet=true;
 }
+
+void 
+ProblemCoreFlows::createKSP()
+{
+       //PETSc Linear solver
+       KSPCreate(PETSC_COMM_WORLD, &_ksp);
+       KSPSetType(_ksp, _ksptype);
+       KSPSetTolerances(_ksp,_precision,_precision,PETSC_DEFAULT,_maxPetscIts);
+       KSPGetPC(_ksp, &_pc);
+       //PETSc preconditioner
+       if(_mpi_size==1 )
+               PCSetType(_pc, _pctype);
+       else
+       {
+               PCSetType(_pc, PCBJACOBI);//Global preconditioner is block jacobi
+               if(_pctype != (char*)&PCILU)//Default pc type is ilu
+               {
+                       PetscOptionsSetValue(NULL,"-sub_pc_type ",_pctype);
+                       PetscOptionsSetValue(NULL,"-sub_ksp_type ","preonly");
+                       //If the above setvalue does not work, try the following
+                       /*
+                       KSPSetUp(_ksp);//to set the block Jacobi data structures (including creation of an internal KSP context for each block)
+                       KSP * subKSP;
+                       PC subpc;
+                       int nlocal;//nb local blocs (should equal 1)
+                       PCBJacobiGetSubKSP(_pc,&nlocal,NULL,&subKSP);
+                       if(nlocal==1)
+                       {
+                               KSPSetType(subKSP[0], KSPPREONLY);//local block solver is same as global
+                               KSPGetPC(subKSP[0],&subpc);
+                               PCSetType(subpc,_pctype);
+                       }
+                       else
+                               throw CdmathException("PC Block Jacobi, more than one block in this processor!!");
+                       */ 
+               }
+       }
+       
+}