From 69717d35b732254023c8fe61648db140b3c1e9f0 Mon Sep 17 00:00:00 2001 From: michael Date: Thu, 13 Jan 2022 18:28:54 +0100 Subject: [PATCH] Added options to visualise matrix --- CoreFlows/Models/inc/ProblemCoreFlows.hxx | 11 +++--- CoreFlows/Models/src/ProblemCoreFlows.cxx | 47 +++++++++++++++++++++-- 2 files changed, 49 insertions(+), 9 deletions(-) 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!!"); + */ + } + } + +} -- 2.39.2