| Directory: | ./ |
|---|---|
| File: | PETScInterface/KSPInterface.cpp |
| Date: | 2024-04-14 07:32:34 |
| Exec | Total | Coverage | |
|---|---|---|---|
| Lines: | 61 | 92 | 66.3% |
| Branches: | 30 | 115 | 26.1% |
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // ______ ______ _ _ _____ ______ | ||
| 2 | // | ____| ____| | (_)/ ____| | ____| | ||
| 3 | // | |__ | |__ | | _| (___ ___| |__ | ||
| 4 | // | __| | __| | | | |\___ \ / __| __| | ||
| 5 | // | | | |____| |____| |____) | (__| |____ | ||
| 6 | // |_| |______|______|_|_____/ \___|______| | ||
| 7 | // Finite Elements for Life Sciences and Engineering | ||
| 8 | // | ||
| 9 | // License: LGL2.1 License | ||
| 10 | // FELiScE default license: LICENSE in root folder | ||
| 11 | // | ||
| 12 | // Main authors: Vicente Mataix Ferrandiz | ||
| 13 | // | ||
| 14 | |||
| 15 | // System includes | ||
| 16 | |||
| 17 | // External includes | ||
| 18 | |||
| 19 | // Project includes | ||
| 20 | #include "Core/mpiInfo.hpp" | ||
| 21 | #include "PETScInterface/KSPInterface.hpp" | ||
| 22 | |||
| 23 | namespace felisce | ||
| 24 | { | ||
| 25 | |||
| 26 | 521 | KSPInterface::KSPInterface() | |
| 27 | { | ||
| 28 | 521 | } | |
| 29 | |||
| 30 | /***********************************************************************************/ | ||
| 31 | /***********************************************************************************/ | ||
| 32 | |||
| 33 | 521 | KSPInterface::~KSPInterface() | |
| 34 | { | ||
| 35 |
2/2✓ Branch 0 taken 230 times.
✓ Branch 1 taken 291 times.
|
521 | if (m_Initialized) { |
| 36 | 230 | KSPDestroy(&m_KSP); | |
| 37 | } | ||
| 38 | 521 | } | |
| 39 | |||
| 40 | /***********************************************************************************/ | ||
| 41 | /***********************************************************************************/ | ||
| 42 | |||
| 43 | 241 | void KSPInterface::init() | |
| 44 | { | ||
| 45 | 241 | init( MpiInfo::petscComm() ); | |
| 46 | 241 | } | |
| 47 | |||
| 48 | /***********************************************************************************/ | ||
| 49 | /***********************************************************************************/ | ||
| 50 | |||
| 51 | 241 | void KSPInterface::init(MPI_Comm comm) | |
| 52 | { | ||
| 53 | 241 | m_comm = comm; | |
| 54 | 241 | KSPCreate(m_comm, &m_KSP); | |
| 55 | 241 | m_Initialized = true; | |
| 56 | 241 | } | |
| 57 | |||
| 58 | /***********************************************************************************/ | ||
| 59 | /***********************************************************************************/ | ||
| 60 | |||
| 61 | 283 | void KSPInterface::initFromSNES(SNESInterface::Pointer snes) | |
| 62 | { | ||
| 63 | 283 | m_comm = snes->m_comm; | |
| 64 | 283 | snes->getKSP(&m_KSP); | |
| 65 | 283 | } | |
| 66 | |||
| 67 | /***********************************************************************************/ | ||
| 68 | /***********************************************************************************/ | ||
| 69 | |||
| 70 | 524 | void KSPInterface::setKSPandPCType( | |
| 71 | const std::string solver, | ||
| 72 | const std::string preconditioner | ||
| 73 | ) | ||
| 74 | { | ||
| 75 |
7/8✓ Branch 1 taken 170 times.
✓ Branch 2 taken 354 times.
✓ Branch 4 taken 170 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4 times.
✓ Branch 8 taken 166 times.
✓ Branch 9 taken 358 times.
✓ Branch 10 taken 166 times.
|
524 | if ( solver == "preonly" || solver == "mumps" || solver == "superlu_dist" ) |
| 76 | 358 | KSPSetType(m_KSP, KSPPREONLY); | |
| 77 | else | ||
| 78 | 166 | KSPSetType(m_KSP, solver.c_str()); | |
| 79 | |||
| 80 | 524 | KSPGetPC(m_KSP,&m_PC); | |
| 81 | 524 | PCSetType(m_PC, preconditioner.c_str()); | |
| 82 | |||
| 83 | 524 | PCSetFromOptions(m_PC); | |
| 84 | 524 | KSPSetFromOptions(m_KSP); | |
| 85 | 524 | } | |
| 86 | |||
| 87 | /***********************************************************************************/ | ||
| 88 | /***********************************************************************************/ | ||
| 89 | |||
| 90 | 524 | void KSPInterface::setTolerances( | |
| 91 | const double relativeTolerance, | ||
| 92 | const double absoluteTolerance, | ||
| 93 | const int maxIteration, | ||
| 94 | const int gmresRestart | ||
| 95 | ) | ||
| 96 | { | ||
| 97 | 524 | KSPSetTolerances(m_KSP,relativeTolerance, absoluteTolerance, PETSC_DEFAULT, maxIteration); | |
| 98 | 524 | KSPGMRESSetRestart(m_KSP, gmresRestart); // ignored by petsc if solver != gmres | |
| 99 | 524 | } | |
| 100 | |||
| 101 | /***********************************************************************************/ | ||
| 102 | /***********************************************************************************/ | ||
| 103 | |||
| 104 | 524 | void KSPInterface::setKSPOptions(const std::string solver, const bool usePrevSol) | |
| 105 | { | ||
| 106 |
5/6✓ Branch 1 taken 170 times.
✓ Branch 2 taken 354 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 170 times.
✓ Branch 6 taken 354 times.
✓ Branch 7 taken 170 times.
|
524 | if ( solver == "preonly" || solver == "mumps" ) { |
| 107 | 354 | PCFactorSetMatSolverType(m_PC, MATSOLVERMUMPS); | |
| 108 | } else | ||
| 109 | #if SUPERLU_DIST_COMPILED | ||
| 110 |
2/2✓ Branch 1 taken 4 times.
✓ Branch 2 taken 166 times.
|
170 | if ( solver == "superlu_dist" ) { |
| 111 | 4 | PCFactorSetMatSolverType(m_PC, MATSOLVERSUPERLU_DIST); | |
| 112 | } else | ||
| 113 | #else | ||
| 114 | if ( solver == "superlu_dist" ) { | ||
| 115 | FEL_WARNING("SuperLUDist not available, we proceed using MUMPS") | ||
| 116 | PCFactorSetMatSolverType(m_PC, MATSOLVERMUMPS); | ||
| 117 | } else | ||
| 118 | #endif | ||
| 119 | { | ||
| 120 |
2/2✓ Branch 0 taken 54 times.
✓ Branch 1 taken 112 times.
|
166 | if ( usePrevSol ) { |
| 121 | 54 | KSPSetInitialGuessNonzero(m_KSP, PETSC_TRUE); | |
| 122 | } | ||
| 123 | } | ||
| 124 | |||
| 125 | { | ||
| 126 | // TODO: to be removed | ||
| 127 | 524 | PetscOptionsSetValue(NULL, "-mat_mumps_icntl_20", "0"); | |
| 128 | } | ||
| 129 | |||
| 130 | 524 | PCSetFromOptions(m_PC); | |
| 131 | 524 | KSPSetFromOptions(m_KSP); | |
| 132 | 524 | } | |
| 133 | |||
| 134 | /***********************************************************************************/ | ||
| 135 | /***********************************************************************************/ | ||
| 136 | |||
| 137 | 6822 | void KSPInterface::setKSPOperator(PetscMatrix& mat, const std::string preconditionerOption) | |
| 138 | { | ||
| 139 | 6822 | PCFactorSetReuseFill(m_PC, PETSC_TRUE); | |
| 140 | |||
| 141 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 6822 times.
|
6822 | if ( preconditionerOption == "SAME_PRECONDITIONER" ) { |
| 142 | ✗ | KSPSetReusePreconditioner(m_KSP, PETSC_TRUE); | |
| 143 | ✗ | PCSetReusePreconditioner(m_PC, PETSC_TRUE); | |
| 144 | } else { | ||
| 145 | 6822 | KSPSetReusePreconditioner(m_KSP, PETSC_FALSE); | |
| 146 | 6822 | PCSetReusePreconditioner(m_PC, PETSC_FALSE); | |
| 147 | } | ||
| 148 | |||
| 149 | 6822 | KSPSetOperators(m_KSP, mat.toPetsc(), mat.toPetsc()); | |
| 150 | 6822 | KSPSetFromOptions(m_KSP); | |
| 151 | 6822 | } | |
| 152 | |||
| 153 | /***********************************************************************************/ | ||
| 154 | /***********************************************************************************/ | ||
| 155 | |||
| 156 | 7806 | void KSPInterface::solve(PetscVector& rRhs, PetscVector& rSol, const int verbose) | |
| 157 | { | ||
| 158 |
3/6✓ Branch 1 taken 7806 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7806 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 7806 times.
✗ Branch 8 not taken.
|
7806 | KSPSolve(m_KSP, rRhs.toPetsc(), rSol.toPetsc()); |
| 159 | |||
| 160 | KSPConvergedReason cvg_reason; | ||
| 161 |
1/2✓ Branch 1 taken 7806 times.
✗ Branch 2 not taken.
|
7806 | KSPGetConvergedReason(m_KSP, &cvg_reason); |
| 162 | PetscInt its; | ||
| 163 |
1/2✓ Branch 1 taken 7806 times.
✗ Branch 2 not taken.
|
7806 | KSPGetIterationNumber(m_KSP, &its); |
| 164 | |||
| 165 |
1/2✓ Branch 0 taken 7806 times.
✗ Branch 1 not taken.
|
7806 | if (cvg_reason > 0) { |
| 166 | |||
| 167 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 7806 times.
|
7806 | if (verbose > 0) { |
| 168 | ✗ | switch (cvg_reason) { | |
| 169 | ✗ | case KSP_CONVERGED_RTOL_NORMAL: PetscPrintf(m_comm, "KSP converged in %d iterations due to ||R|| < rtol*||R_initial|| (NORMAL)\n", its); break; | |
| 170 | ✗ | case KSP_CONVERGED_ATOL_NORMAL: PetscPrintf(m_comm, "KSP converged in %d iterations due to ||R|| < atol (NORMAL)\n", its); break; | |
| 171 | ✗ | case KSP_CONVERGED_RTOL: PetscPrintf(m_comm, "KSP converged in %d iterations due to ||R|| < rtol*||R_initial|| \n", its); break; | |
| 172 | ✗ | case KSP_CONVERGED_ATOL: PetscPrintf(m_comm, "KSP converged in %d iterations due to ||R|| < atol \n", its); break; | |
| 173 | ✗ | case KSP_CONVERGED_ITS: PetscPrintf(m_comm, "KSP converged in %d iteration, direct solver \n", its); break; | |
| 174 | #if PETSC_VERSION_LT(3, 19, 0) | ||
| 175 | case KSP_CONVERGED_CG_NEG_CURVE: PetscPrintf(m_comm, "KSP converged in %d iterations, due to KSP_CONVERGED_CG_NEG_CURVE, returned only by KSPNASH used in SNES solver \n", its); break; | ||
| 176 | case KSP_CONVERGED_CG_CONSTRAINED: PetscPrintf(m_comm, "KSP converged in %d iterations, due to KSP_CONVERGED_CG_CONSTRAINED, returned only by KSPSTCG used in SNES solver \n", its); break; | ||
| 177 | case KSP_CONVERGED_STEP_LENGTH: PetscPrintf(m_comm, "KSP converged in %d iterations, due to KSP_CONVERGED_STEP_LENGTH, returned only by KSPLGTR used in SNES solver \n", its); break; | ||
| 178 | #else | ||
| 179 | ✗ | case KSP_CONVERGED_NEG_CURVE: PetscPrintf(m_comm, "KSP converged in %d iterations, due to KSP_CONVERGED_NEG_CURVE, returned only by KSPNASH, KSPSTCG and KSPLGTR used in SNES solver \n", its); break; | |
| 180 | ✗ | case KSP_CONVERGED_STEP_LENGTH: PetscPrintf(m_comm, "KSP converged in %d iterations, due to KSP_CONVERGED_STEP_LENGTH, returned only by KSPNASH, KSPSTCG and KSPLGTR used in SNES solver \n", its); break; | |
| 181 | #endif | ||
| 182 | ✗ | case KSP_CONVERGED_HAPPY_BREAKDOWN: PetscPrintf(m_comm, "KSP converged in %d iterations due to happy breakdown \n", its); break; | |
| 183 | ✗ | default: PetscPrintf(m_comm, "KSP converged in %d iterations \n", its); | |
| 184 | } | ||
| 185 | } | ||
| 186 | |||
| 187 | KSPType ksptype; | ||
| 188 |
1/2✓ Branch 1 taken 7806 times.
✗ Branch 2 not taken.
|
7806 | KSPGetType(m_KSP, &ksptype); |
| 189 | PCType pctype; | ||
| 190 |
1/2✓ Branch 1 taken 7806 times.
✗ Branch 2 not taken.
|
7806 | PCGetType(m_PC, &pctype); |
| 191 | PetscReal rnorm; | ||
| 192 |
1/2✓ Branch 1 taken 7806 times.
✗ Branch 2 not taken.
|
7806 | KSPGetResidualNorm(m_KSP, &rnorm); |
| 193 | |||
| 194 |
1/2✓ Branch 1 taken 7806 times.
✗ Branch 2 not taken.
|
7806 | PetscPrintf(m_comm, "KSP %s-%s converged in %d iterations (residual = %e)\n", ksptype, pctype, its, rnorm); |
| 195 | } else { | ||
| 196 | ✗ | std::string error_message = "Convergence failed in KSP solver. "; | |
| 197 | ✗ | switch (cvg_reason) { | |
| 198 | ✗ | case KSP_CONVERGED_ITERATING: PetscPrintf(m_comm, "KSP failed due to a call to KSPGetConvergedReason() while KSPSolve() is still running \n"); break; | |
| 199 | ✗ | case KSP_DIVERGED_NULL: PetscPrintf(m_comm, "KSP failed due to KSP_DIVERGED_NULL \n"); break; | |
| 200 | ✗ | case KSP_DIVERGED_ITS: PetscPrintf(m_comm, "KSP diverged due to KSPSolve() has reached the maximum number of iterations %d \n", its); break; | |
| 201 | ✗ | case KSP_DIVERGED_DTOL: PetscPrintf(m_comm, "KSP diverged due to ||R|| > divtol*||R_initial|| \n"); break; | |
| 202 | ✗ | case KSP_DIVERGED_BREAKDOWN: PetscPrintf(m_comm, "KSP failed due to a generic breakdown in method. See KSP_DIVERGED_BREAKDOWN reason \n"); break; | |
| 203 | ✗ | case KSP_DIVERGED_BREAKDOWN_BICG: PetscPrintf(m_comm, "KSP failed due to a generic breakdown in BICG method. See KSP_DIVERGED_BREAKDOWN_BICG reason \n"); break; | |
| 204 | ✗ | case KSP_DIVERGED_NONSYMMETRIC: PetscPrintf(m_comm, "KSP failed due to matrix or preconditioner not symmetric while the Krylov method requires symmetry \n"); break; | |
| 205 | ✗ | case KSP_DIVERGED_INDEFINITE_PC: PetscPrintf(m_comm, "KSP failed due to indefinite preconditioner (with both positive and negative eigenvalues) while the Krylov method requires it be positive definite \n"); break; | |
| 206 | ✗ | case KSP_DIVERGED_NANORINF: PetscPrintf(m_comm, "KSP failed due to residual norm equal to Not-a-number or Inf \n"); break; | |
| 207 | ✗ | case KSP_DIVERGED_INDEFINITE_MAT: PetscPrintf(m_comm, "KSP failed due to indefinite matrix (with both positive and negative eigenvalues) while the Krylov method requires it be positive definite \n"); break; | |
| 208 | ✗ | case KSP_DIVERGED_PC_FAILED: PetscPrintf(m_comm, "KSP failed due to was not possible to build or use the requested preconditioner (usually due to zero pivot in a factorization) \n"); break; | |
| 209 | ✗ | default: PetscPrintf(m_comm, "\n"); | |
| 210 | } | ||
| 211 | ✗ | FEL_ERROR(error_message); | |
| 212 | } | ||
| 213 | 7806 | } | |
| 214 | |||
| 215 | /***********************************************************************************/ | ||
| 216 | /***********************************************************************************/ | ||
| 217 | |||
| 218 | ✗ | void KSPInterface::destroy() | |
| 219 | { | ||
| 220 | ✗ | if (m_Initialized) { | |
| 221 | ✗ | KSPDestroy(&m_KSP); | |
| 222 | ✗ | m_Initialized = false; | |
| 223 | } | ||
| 224 | } | ||
| 225 | |||
| 226 | } | ||
| 227 |