From: michael <michael@localhost.localdomain>
Date: Thu, 13 Jan 2022 17:28:54 +0000 (+0100)
Subject: Added options to visualise matrix
X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=69717d35b732254023c8fe61648db140b3c1e9f0;p=tools%2Fsolverlab.git

Added options to visualise matrix
---

diff --git a/CoreFlows/Models/inc/ProblemCoreFlows.hxx b/CoreFlows/Models/inc/ProblemCoreFlows.hxx
index ee2245e..738f049 100755
--- 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 6c35c63..492dc34 100755
--- 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!!");
+			*/ 
+		}
+	}
+	
+}