1 /* 2 Defines the multigrid preconditioner interface. 3 */ 4 #include <petsc/private/pcmgimpl.h> /*I "petscksp.h" I*/ 5 #include <petsc/private/kspimpl.h> 6 #include <petscdm.h> 7 PETSC_INTERN PetscErrorCode PCPreSolveChangeRHS(PC, PetscBool *); 8 9 /* 10 Contains the list of registered coarse space construction routines 11 */ 12 PetscFunctionList PCMGCoarseList = NULL; 13 14 PetscErrorCode PCMGMCycle_Private(PC pc, PC_MG_Levels **mglevelsin, PetscBool transpose, PetscBool matapp, PCRichardsonConvergedReason *reason) 15 { 16 PC_MG *mg = (PC_MG *)pc->data; 17 PC_MG_Levels *mgc, *mglevels = *mglevelsin; 18 PetscInt cycles = (mglevels->level == 1) ? 1 : mglevels->cycles; 19 20 PetscFunctionBegin; 21 if (mglevels->eventsmoothsolve) PetscCall(PetscLogEventBegin(mglevels->eventsmoothsolve, 0, 0, 0, 0)); 22 if (!transpose) { 23 if (matapp) { 24 PetscCall(KSPMatSolve(mglevels->smoothd, mglevels->B, mglevels->X)); /* pre-smooth */ 25 PetscCall(KSPCheckSolve(mglevels->smoothd, pc, NULL)); 26 } else { 27 PetscCall(KSPSolve(mglevels->smoothd, mglevels->b, mglevels->x)); /* pre-smooth */ 28 PetscCall(KSPCheckSolve(mglevels->smoothd, pc, mglevels->x)); 29 } 30 } else { 31 PetscCheck(!matapp, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Not supported"); 32 PetscCall(KSPSolveTranspose(mglevels->smoothu, mglevels->b, mglevels->x)); /* transpose of post-smooth */ 33 PetscCall(KSPCheckSolve(mglevels->smoothu, pc, mglevels->x)); 34 } 35 if (mglevels->eventsmoothsolve) PetscCall(PetscLogEventEnd(mglevels->eventsmoothsolve, 0, 0, 0, 0)); 36 if (mglevels->level) { /* not the coarsest grid */ 37 if (mglevels->eventresidual) PetscCall(PetscLogEventBegin(mglevels->eventresidual, 0, 0, 0, 0)); 38 if (matapp && !mglevels->R) PetscCall(MatDuplicate(mglevels->B, MAT_DO_NOT_COPY_VALUES, &mglevels->R)); 39 if (!transpose) { 40 if (matapp) PetscCall((*mglevels->matresidual)(mglevels->A, mglevels->B, mglevels->X, mglevels->R)); 41 else PetscCall((*mglevels->residual)(mglevels->A, mglevels->b, mglevels->x, mglevels->r)); 42 } else { 43 if (matapp) PetscCall((*mglevels->matresidualtranspose)(mglevels->A, mglevels->B, mglevels->X, mglevels->R)); 44 else PetscCall((*mglevels->residualtranspose)(mglevels->A, mglevels->b, mglevels->x, mglevels->r)); 45 } 46 if (mglevels->eventresidual) PetscCall(PetscLogEventEnd(mglevels->eventresidual, 0, 0, 0, 0)); 47 48 /* if on finest level and have convergence criteria set */ 49 if (mglevels->level == mglevels->levels - 1 && mg->ttol && reason) { 50 PetscReal rnorm; 51 PetscCall(VecNorm(mglevels->r, NORM_2, &rnorm)); 52 if (rnorm <= mg->ttol) { 53 if (rnorm < mg->abstol) { 54 *reason = PCRICHARDSON_CONVERGED_ATOL; 55 PetscCall(PetscInfo(pc, "Linear solver has converged. Residual norm %g is less than absolute tolerance %g\n", (double)rnorm, (double)mg->abstol)); 56 } else { 57 *reason = PCRICHARDSON_CONVERGED_RTOL; 58 PetscCall(PetscInfo(pc, "Linear solver has converged. Residual norm %g is less than relative tolerance times initial residual norm %g\n", (double)rnorm, (double)mg->ttol)); 59 } 60 PetscFunctionReturn(PETSC_SUCCESS); 61 } 62 } 63 64 mgc = *(mglevelsin - 1); 65 if (mglevels->eventinterprestrict) PetscCall(PetscLogEventBegin(mglevels->eventinterprestrict, 0, 0, 0, 0)); 66 if (!transpose) { 67 if (matapp) PetscCall(MatMatRestrict(mglevels->restrct, mglevels->R, &mgc->B)); 68 else PetscCall(MatRestrict(mglevels->restrct, mglevels->r, mgc->b)); 69 } else { 70 if (matapp) PetscCall(MatMatRestrict(mglevels->interpolate, mglevels->R, &mgc->B)); 71 else PetscCall(MatRestrict(mglevels->interpolate, mglevels->r, mgc->b)); 72 } 73 if (mglevels->eventinterprestrict) PetscCall(PetscLogEventEnd(mglevels->eventinterprestrict, 0, 0, 0, 0)); 74 if (matapp) { 75 if (!mgc->X) { 76 PetscCall(MatDuplicate(mgc->B, MAT_DO_NOT_COPY_VALUES, &mgc->X)); 77 } else { 78 PetscCall(MatZeroEntries(mgc->X)); 79 } 80 } else { 81 PetscCall(VecZeroEntries(mgc->x)); 82 } 83 while (cycles--) PetscCall(PCMGMCycle_Private(pc, mglevelsin - 1, transpose, matapp, reason)); 84 if (mglevels->eventinterprestrict) PetscCall(PetscLogEventBegin(mglevels->eventinterprestrict, 0, 0, 0, 0)); 85 if (!transpose) { 86 if (matapp) PetscCall(MatMatInterpolateAdd(mglevels->interpolate, mgc->X, mglevels->X, &mglevels->X)); 87 else PetscCall(MatInterpolateAdd(mglevels->interpolate, mgc->x, mglevels->x, mglevels->x)); 88 } else { 89 PetscCall(MatInterpolateAdd(mglevels->restrct, mgc->x, mglevels->x, mglevels->x)); 90 } 91 if (mglevels->eventinterprestrict) PetscCall(PetscLogEventEnd(mglevels->eventinterprestrict, 0, 0, 0, 0)); 92 if (mglevels->eventsmoothsolve) PetscCall(PetscLogEventBegin(mglevels->eventsmoothsolve, 0, 0, 0, 0)); 93 if (!transpose) { 94 if (matapp) { 95 PetscCall(KSPMatSolve(mglevels->smoothu, mglevels->B, mglevels->X)); /* post smooth */ 96 PetscCall(KSPCheckSolve(mglevels->smoothu, pc, NULL)); 97 } else { 98 PetscCall(KSPSolve(mglevels->smoothu, mglevels->b, mglevels->x)); /* post smooth */ 99 PetscCall(KSPCheckSolve(mglevels->smoothu, pc, mglevels->x)); 100 } 101 } else { 102 PetscCheck(!matapp, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Not supported"); 103 PetscCall(KSPSolveTranspose(mglevels->smoothd, mglevels->b, mglevels->x)); /* post smooth */ 104 PetscCall(KSPCheckSolve(mglevels->smoothd, pc, mglevels->x)); 105 } 106 if (mglevels->cr) { 107 PetscCheck(!matapp, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Not supported"); 108 /* TODO Turn on copy and turn off noisy if we have an exact solution 109 PetscCall(VecCopy(mglevels->x, mglevels->crx)); 110 PetscCall(VecCopy(mglevels->b, mglevels->crb)); */ 111 PetscCall(KSPSetNoisy_Private(mglevels->crx)); 112 PetscCall(KSPSolve(mglevels->cr, mglevels->crb, mglevels->crx)); /* compatible relaxation */ 113 PetscCall(KSPCheckSolve(mglevels->cr, pc, mglevels->crx)); 114 } 115 if (mglevels->eventsmoothsolve) PetscCall(PetscLogEventEnd(mglevels->eventsmoothsolve, 0, 0, 0, 0)); 116 } 117 PetscFunctionReturn(PETSC_SUCCESS); 118 } 119 120 static PetscErrorCode PCApplyRichardson_MG(PC pc, Vec b, Vec x, Vec w, PetscReal rtol, PetscReal abstol, PetscReal dtol, PetscInt its, PetscBool zeroguess, PetscInt *outits, PCRichardsonConvergedReason *reason) 121 { 122 PC_MG *mg = (PC_MG *)pc->data; 123 PC_MG_Levels **mglevels = mg->levels; 124 PC tpc; 125 PetscBool changeu, changed; 126 PetscInt levels = mglevels[0]->levels, i; 127 128 PetscFunctionBegin; 129 /* When the DM is supplying the matrix then it will not exist until here */ 130 for (i = 0; i < levels; i++) { 131 if (!mglevels[i]->A) { 132 PetscCall(KSPGetOperators(mglevels[i]->smoothu, &mglevels[i]->A, NULL)); 133 PetscCall(PetscObjectReference((PetscObject)mglevels[i]->A)); 134 } 135 } 136 137 PetscCall(KSPGetPC(mglevels[levels - 1]->smoothd, &tpc)); 138 PetscCall(PCPreSolveChangeRHS(tpc, &changed)); 139 PetscCall(KSPGetPC(mglevels[levels - 1]->smoothu, &tpc)); 140 PetscCall(PCPreSolveChangeRHS(tpc, &changeu)); 141 if (!changed && !changeu) { 142 PetscCall(VecDestroy(&mglevels[levels - 1]->b)); 143 mglevels[levels - 1]->b = b; 144 } else { /* if the smoother changes the rhs during PreSolve, we cannot use the input vector */ 145 if (!mglevels[levels - 1]->b) { 146 Vec *vec; 147 148 PetscCall(KSPCreateVecs(mglevels[levels - 1]->smoothd, 1, &vec, 0, NULL)); 149 mglevels[levels - 1]->b = *vec; 150 PetscCall(PetscFree(vec)); 151 } 152 PetscCall(VecCopy(b, mglevels[levels - 1]->b)); 153 } 154 mglevels[levels - 1]->x = x; 155 156 mg->rtol = rtol; 157 mg->abstol = abstol; 158 mg->dtol = dtol; 159 if (rtol) { 160 /* compute initial residual norm for relative convergence test */ 161 PetscReal rnorm; 162 if (zeroguess) { 163 PetscCall(VecNorm(b, NORM_2, &rnorm)); 164 } else { 165 PetscCall((*mglevels[levels - 1]->residual)(mglevels[levels - 1]->A, b, x, w)); 166 PetscCall(VecNorm(w, NORM_2, &rnorm)); 167 } 168 mg->ttol = PetscMax(rtol * rnorm, abstol); 169 } else if (abstol) mg->ttol = abstol; 170 else mg->ttol = 0.0; 171 172 /* since smoother is applied to full system, not just residual we need to make sure that smoothers don't 173 stop prematurely due to small residual */ 174 for (i = 1; i < levels; i++) { 175 PetscCall(KSPSetTolerances(mglevels[i]->smoothu, 0, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT)); 176 if (mglevels[i]->smoothu != mglevels[i]->smoothd) { 177 /* For Richardson the initial guess is nonzero since it is solving in each cycle the original system not just applying as a preconditioner */ 178 PetscCall(KSPSetInitialGuessNonzero(mglevels[i]->smoothd, PETSC_TRUE)); 179 PetscCall(KSPSetTolerances(mglevels[i]->smoothd, 0, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT)); 180 } 181 } 182 183 *reason = PCRICHARDSON_NOT_SET; 184 for (i = 0; i < its; i++) { 185 PetscCall(PCMGMCycle_Private(pc, mglevels + levels - 1, PETSC_FALSE, PETSC_FALSE, reason)); 186 if (*reason) break; 187 } 188 if (*reason == PCRICHARDSON_NOT_SET) *reason = PCRICHARDSON_CONVERGED_ITS; 189 *outits = i; 190 if (!changed && !changeu) mglevels[levels - 1]->b = NULL; 191 PetscFunctionReturn(PETSC_SUCCESS); 192 } 193 194 PetscErrorCode PCReset_MG(PC pc) 195 { 196 PC_MG *mg = (PC_MG *)pc->data; 197 PC_MG_Levels **mglevels = mg->levels; 198 PetscInt i, n; 199 200 PetscFunctionBegin; 201 if (mglevels) { 202 n = mglevels[0]->levels; 203 for (i = 0; i < n - 1; i++) { 204 PetscCall(VecDestroy(&mglevels[i + 1]->r)); 205 PetscCall(VecDestroy(&mglevels[i]->b)); 206 PetscCall(VecDestroy(&mglevels[i]->x)); 207 PetscCall(MatDestroy(&mglevels[i + 1]->R)); 208 PetscCall(MatDestroy(&mglevels[i]->B)); 209 PetscCall(MatDestroy(&mglevels[i]->X)); 210 PetscCall(VecDestroy(&mglevels[i]->crx)); 211 PetscCall(VecDestroy(&mglevels[i]->crb)); 212 PetscCall(MatDestroy(&mglevels[i + 1]->restrct)); 213 PetscCall(MatDestroy(&mglevels[i + 1]->interpolate)); 214 PetscCall(MatDestroy(&mglevels[i + 1]->inject)); 215 PetscCall(VecDestroy(&mglevels[i + 1]->rscale)); 216 } 217 PetscCall(VecDestroy(&mglevels[n - 1]->crx)); 218 PetscCall(VecDestroy(&mglevels[n - 1]->crb)); 219 /* this is not null only if the smoother on the finest level 220 changes the rhs during PreSolve */ 221 PetscCall(VecDestroy(&mglevels[n - 1]->b)); 222 PetscCall(MatDestroy(&mglevels[n - 1]->B)); 223 224 for (i = 0; i < n; i++) { 225 PetscCall(MatDestroy(&mglevels[i]->coarseSpace)); 226 PetscCall(MatDestroy(&mglevels[i]->A)); 227 if (mglevels[i]->smoothd != mglevels[i]->smoothu) PetscCall(KSPReset(mglevels[i]->smoothd)); 228 PetscCall(KSPReset(mglevels[i]->smoothu)); 229 if (mglevels[i]->cr) PetscCall(KSPReset(mglevels[i]->cr)); 230 } 231 mg->Nc = 0; 232 } 233 PetscFunctionReturn(PETSC_SUCCESS); 234 } 235 236 /* Implementing CR 237 238 We only want to make corrections that ``do not change'' the coarse solution. What we mean by not changing is that if I prolong my coarse solution to the fine grid and then inject that fine solution back to the coarse grid, I get the same answer. Injection is what Brannick calls R. We want the complementary projector to Inj, which we will call S, after Brannick, so that Inj S = 0. Now the orthogonal projector onto the range of Inj^T is 239 240 Inj^T (Inj Inj^T)^{-1} Inj 241 242 and if Inj is a VecScatter, as it is now in PETSc, we have 243 244 Inj^T Inj 245 246 and 247 248 S = I - Inj^T Inj 249 250 since 251 252 Inj S = Inj - (Inj Inj^T) Inj = 0. 253 254 Brannick suggests 255 256 A \to S^T A S \qquad\mathrm{and}\qquad M \to S^T M S 257 258 but I do not think his :math:`S^T S = I` is correct. Our S is an orthogonal projector, so :math:`S^T S = S^2 = S`. We will use 259 260 M^{-1} A \to S M^{-1} A S 261 262 In fact, since it is somewhat hard in PETSc to do the symmetric application, we will just apply S on the left. 263 264 Check: || Inj P - I ||_F < tol 265 Check: In general, Inj Inj^T = I 266 */ 267 268 typedef struct { 269 PC mg; /* The PCMG object */ 270 PetscInt l; /* The multigrid level for this solver */ 271 Mat Inj; /* The injection matrix */ 272 Mat S; /* I - Inj^T Inj */ 273 } CRContext; 274 275 static PetscErrorCode CRSetup_Private(PC pc) 276 { 277 CRContext *ctx; 278 Mat It; 279 280 PetscFunctionBeginUser; 281 PetscCall(PCShellGetContext(pc, &ctx)); 282 PetscCall(PCMGGetInjection(ctx->mg, ctx->l, &It)); 283 PetscCheck(It, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "CR requires that injection be defined for this PCMG"); 284 PetscCall(MatCreateTranspose(It, &ctx->Inj)); 285 PetscCall(MatCreateNormal(ctx->Inj, &ctx->S)); 286 PetscCall(MatScale(ctx->S, -1.0)); 287 PetscCall(MatShift(ctx->S, 1.0)); 288 PetscFunctionReturn(PETSC_SUCCESS); 289 } 290 291 static PetscErrorCode CRApply_Private(PC pc, Vec x, Vec y) 292 { 293 CRContext *ctx; 294 295 PetscFunctionBeginUser; 296 PetscCall(PCShellGetContext(pc, &ctx)); 297 PetscCall(MatMult(ctx->S, x, y)); 298 PetscFunctionReturn(PETSC_SUCCESS); 299 } 300 301 static PetscErrorCode CRDestroy_Private(PC pc) 302 { 303 CRContext *ctx; 304 305 PetscFunctionBeginUser; 306 PetscCall(PCShellGetContext(pc, &ctx)); 307 PetscCall(MatDestroy(&ctx->Inj)); 308 PetscCall(MatDestroy(&ctx->S)); 309 PetscCall(PetscFree(ctx)); 310 PetscCall(PCShellSetContext(pc, NULL)); 311 PetscFunctionReturn(PETSC_SUCCESS); 312 } 313 314 static PetscErrorCode CreateCR_Private(PC pc, PetscInt l, PC *cr) 315 { 316 CRContext *ctx; 317 318 PetscFunctionBeginUser; 319 PetscCall(PCCreate(PetscObjectComm((PetscObject)pc), cr)); 320 PetscCall(PetscObjectSetName((PetscObject)*cr, "S (complementary projector to injection)")); 321 PetscCall(PetscCalloc1(1, &ctx)); 322 ctx->mg = pc; 323 ctx->l = l; 324 PetscCall(PCSetType(*cr, PCSHELL)); 325 PetscCall(PCShellSetContext(*cr, ctx)); 326 PetscCall(PCShellSetApply(*cr, CRApply_Private)); 327 PetscCall(PCShellSetSetUp(*cr, CRSetup_Private)); 328 PetscCall(PCShellSetDestroy(*cr, CRDestroy_Private)); 329 PetscFunctionReturn(PETSC_SUCCESS); 330 } 331 332 PETSC_EXTERN PetscErrorCode PetscOptionsFindPairPrefix_Private(PetscOptions, const char[], const char[], const char *[], const char *[], PetscBool *); 333 334 PetscErrorCode PCMGSetLevels_MG(PC pc, PetscInt levels, MPI_Comm *comms) 335 { 336 PC_MG *mg = (PC_MG *)pc->data; 337 MPI_Comm comm; 338 PC_MG_Levels **mglevels = mg->levels; 339 PCMGType mgtype = mg->am; 340 PetscInt mgctype = (PetscInt)PC_MG_CYCLE_V; 341 PetscInt i; 342 PetscMPIInt size; 343 const char *prefix; 344 PC ipc; 345 PetscInt n; 346 347 PetscFunctionBegin; 348 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 349 PetscValidLogicalCollectiveInt(pc, levels, 2); 350 if (mg->nlevels == levels) PetscFunctionReturn(PETSC_SUCCESS); 351 PetscCall(PetscObjectGetComm((PetscObject)pc, &comm)); 352 if (mglevels) { 353 mgctype = mglevels[0]->cycles; 354 /* changing the number of levels so free up the previous stuff */ 355 PetscCall(PCReset_MG(pc)); 356 n = mglevels[0]->levels; 357 for (i = 0; i < n; i++) { 358 if (mglevels[i]->smoothd != mglevels[i]->smoothu) PetscCall(KSPDestroy(&mglevels[i]->smoothd)); 359 PetscCall(KSPDestroy(&mglevels[i]->smoothu)); 360 PetscCall(KSPDestroy(&mglevels[i]->cr)); 361 PetscCall(PetscFree(mglevels[i])); 362 } 363 PetscCall(PetscFree(mg->levels)); 364 } 365 366 mg->nlevels = levels; 367 368 PetscCall(PetscMalloc1(levels, &mglevels)); 369 370 PetscCall(PCGetOptionsPrefix(pc, &prefix)); 371 372 mg->stageApply = 0; 373 for (i = 0; i < levels; i++) { 374 PetscCall(PetscNew(&mglevels[i])); 375 376 mglevels[i]->level = i; 377 mglevels[i]->levels = levels; 378 mglevels[i]->cycles = mgctype; 379 mg->default_smoothu = 2; 380 mg->default_smoothd = 2; 381 mglevels[i]->eventsmoothsetup = 0; 382 mglevels[i]->eventsmoothsolve = 0; 383 mglevels[i]->eventresidual = 0; 384 mglevels[i]->eventinterprestrict = 0; 385 386 if (comms) comm = comms[i]; 387 if (comm != MPI_COMM_NULL) { 388 PetscCall(KSPCreate(comm, &mglevels[i]->smoothd)); 389 PetscCall(KSPSetNestLevel(mglevels[i]->smoothd, pc->kspnestlevel)); 390 PetscCall(KSPSetErrorIfNotConverged(mglevels[i]->smoothd, pc->erroriffailure)); 391 PetscCall(PetscObjectIncrementTabLevel((PetscObject)mglevels[i]->smoothd, (PetscObject)pc, levels - i)); 392 PetscCall(KSPSetOptionsPrefix(mglevels[i]->smoothd, prefix)); 393 PetscCall(PetscObjectComposedDataSetInt((PetscObject)mglevels[i]->smoothd, PetscMGLevelId, mglevels[i]->level)); 394 if (i == 0 && levels > 1) { // coarse grid 395 PetscCall(KSPAppendOptionsPrefix(mglevels[0]->smoothd, "mg_coarse_")); 396 397 /* coarse solve is (redundant) LU by default; set shifttype NONZERO to avoid annoying zero-pivot in LU preconditioner */ 398 PetscCall(KSPSetType(mglevels[0]->smoothd, KSPPREONLY)); 399 PetscCall(KSPGetPC(mglevels[0]->smoothd, &ipc)); 400 PetscCallMPI(MPI_Comm_size(comm, &size)); 401 if (size > 1) { 402 PetscCall(PCSetType(ipc, PCREDUNDANT)); 403 } else { 404 PetscCall(PCSetType(ipc, PCLU)); 405 } 406 PetscCall(PCFactorSetShiftType(ipc, MAT_SHIFT_INBLOCKS)); 407 } else { 408 char tprefix[128]; 409 410 PetscCall(KSPSetType(mglevels[i]->smoothd, KSPCHEBYSHEV)); 411 PetscCall(KSPSetConvergenceTest(mglevels[i]->smoothd, KSPConvergedSkip, NULL, NULL)); 412 PetscCall(KSPSetNormType(mglevels[i]->smoothd, KSP_NORM_NONE)); 413 PetscCall(KSPGetPC(mglevels[i]->smoothd, &ipc)); 414 PetscCall(PCSetType(ipc, PCSOR)); 415 PetscCall(KSPSetTolerances(mglevels[i]->smoothd, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT, mg->default_smoothd)); 416 417 if (i == levels - 1 && levels > 1) { // replace 'mg_finegrid_' with 'mg_levels_X_' 418 PetscBool set; 419 PetscCall(PetscOptionsFindPairPrefix_Private(((PetscObject)mglevels[i]->smoothd)->options, ((PetscObject)mglevels[i]->smoothd)->prefix, "-mg_fine_", NULL, NULL, &set)); 420 if (set) { 421 if (prefix) PetscCall(PetscSNPrintf(tprefix, 128, "%smg_fine_", prefix)); 422 else PetscCall(PetscSNPrintf(tprefix, 128, "mg_fine_")); 423 PetscCall(KSPSetOptionsPrefix(mglevels[i]->smoothd, tprefix)); 424 } else { 425 PetscCall(PetscSNPrintf(tprefix, 128, "mg_levels_%" PetscInt_FMT "_", i)); 426 PetscCall(KSPAppendOptionsPrefix(mglevels[i]->smoothd, tprefix)); 427 } 428 } else { 429 PetscCall(PetscSNPrintf(tprefix, 128, "mg_levels_%" PetscInt_FMT "_", i)); 430 PetscCall(KSPAppendOptionsPrefix(mglevels[i]->smoothd, tprefix)); 431 } 432 } 433 } 434 mglevels[i]->smoothu = mglevels[i]->smoothd; 435 mg->rtol = 0.0; 436 mg->abstol = 0.0; 437 mg->dtol = 0.0; 438 mg->ttol = 0.0; 439 mg->cyclesperpcapply = 1; 440 } 441 mg->levels = mglevels; 442 PetscCall(PCMGSetType(pc, mgtype)); 443 PetscFunctionReturn(PETSC_SUCCESS); 444 } 445 446 /*@C 447 PCMGSetLevels - Sets the number of levels to use with `PCMG`. 448 Must be called before any other `PCMG` routine. 449 450 Logically Collective 451 452 Input Parameters: 453 + pc - the preconditioner context 454 . levels - the number of levels 455 - comms - optional communicators for each level; this is to allow solving the coarser problems 456 on smaller sets of processes. For processes that are not included in the computation 457 you must pass `MPI_COMM_NULL`. Use comms = `NULL` to specify that all processes 458 should participate in each level of problem. 459 460 Level: intermediate 461 462 Notes: 463 If the number of levels is one then the multigrid uses the `-mg_levels` prefix 464 for setting the level options rather than the `-mg_coarse` or `-mg_fine` prefix. 465 466 You can free the information in comms after this routine is called. 467 468 The array of MPI communicators must contain `MPI_COMM_NULL` for those ranks that at each level 469 are not participating in the coarser solve. For example, with 2 levels and 1 and 2 ranks on 470 the two levels, rank 0 in the original communicator will pass in an array of 2 communicators 471 of size 2 and 1, while rank 1 in the original communicator will pass in array of 2 communicators 472 the first of size 2 and the second of value `MPI_COMM_NULL` since the rank 1 does not participate 473 in the coarse grid solve. 474 475 Since each coarser level may have a new `MPI_Comm` with fewer ranks than the previous, one 476 must take special care in providing the restriction and interpolation operation. We recommend 477 providing these as two step operations; first perform a standard restriction or interpolation on 478 the full number of ranks for that level and then use an MPI call to copy the resulting vector 479 array entries (after calls to VecGetArray()) to the smaller or larger number of ranks, note in both 480 cases the MPI calls must be made on the larger of the two communicators. Traditional MPI send and 481 receives or `MPI_AlltoAllv()` could be used to do the reshuffling of the vector entries. 482 483 Fortran Notes: 484 Use comms = `PETSC_NULL_MPI_COMM` as the equivalent of `NULL` in the C interface. Note `PETSC_NULL_MPI_COMM` 485 is not `MPI_COMM_NULL`. It is more like `PETSC_NULL_INTEGER`, `PETSC_NULL_REAL` etc. 486 487 .seealso: [](ch_ksp), `PCMGSetType()`, `PCMGGetLevels()` 488 @*/ 489 PetscErrorCode PCMGSetLevels(PC pc, PetscInt levels, MPI_Comm *comms) 490 { 491 PetscFunctionBegin; 492 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 493 if (comms) PetscAssertPointer(comms, 3); 494 PetscTryMethod(pc, "PCMGSetLevels_C", (PC, PetscInt, MPI_Comm *), (pc, levels, comms)); 495 PetscFunctionReturn(PETSC_SUCCESS); 496 } 497 498 PetscErrorCode PCDestroy_MG(PC pc) 499 { 500 PC_MG *mg = (PC_MG *)pc->data; 501 PC_MG_Levels **mglevels = mg->levels; 502 PetscInt i, n; 503 504 PetscFunctionBegin; 505 PetscCall(PCReset_MG(pc)); 506 if (mglevels) { 507 n = mglevels[0]->levels; 508 for (i = 0; i < n; i++) { 509 if (mglevels[i]->smoothd != mglevels[i]->smoothu) PetscCall(KSPDestroy(&mglevels[i]->smoothd)); 510 PetscCall(KSPDestroy(&mglevels[i]->smoothu)); 511 PetscCall(KSPDestroy(&mglevels[i]->cr)); 512 PetscCall(PetscFree(mglevels[i])); 513 } 514 PetscCall(PetscFree(mg->levels)); 515 } 516 PetscCall(PetscFree(pc->data)); 517 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGetInterpolations_C", NULL)); 518 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGetCoarseOperators_C", NULL)); 519 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGSetGalerkin_C", NULL)); 520 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGetLevels_C", NULL)); 521 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGSetLevels_C", NULL)); 522 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGetInterpolations_C", NULL)); 523 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGetCoarseOperators_C", NULL)); 524 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGSetAdaptInterpolation_C", NULL)); 525 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGetAdaptInterpolation_C", NULL)); 526 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGSetAdaptCR_C", NULL)); 527 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGetAdaptCR_C", NULL)); 528 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGSetAdaptCoarseSpaceType_C", NULL)); 529 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGetAdaptCoarseSpaceType_C", NULL)); 530 PetscFunctionReturn(PETSC_SUCCESS); 531 } 532 533 /* 534 PCApply_MG - Runs either an additive, multiplicative, Kaskadic 535 or full cycle of multigrid. 536 537 Note: 538 A simple wrapper which calls PCMGMCycle(),PCMGACycle(), or PCMGFCycle(). 539 */ 540 static PetscErrorCode PCApply_MG_Internal(PC pc, Vec b, Vec x, Mat B, Mat X, PetscBool transpose) 541 { 542 PC_MG *mg = (PC_MG *)pc->data; 543 PC_MG_Levels **mglevels = mg->levels; 544 PC tpc; 545 PetscInt levels = mglevels[0]->levels, i; 546 PetscBool changeu, changed, matapp; 547 548 PetscFunctionBegin; 549 matapp = (PetscBool)(B && X); 550 if (mg->stageApply) PetscCall(PetscLogStagePush(mg->stageApply)); 551 /* When the DM is supplying the matrix then it will not exist until here */ 552 for (i = 0; i < levels; i++) { 553 if (!mglevels[i]->A) { 554 PetscCall(KSPGetOperators(mglevels[i]->smoothu, &mglevels[i]->A, NULL)); 555 PetscCall(PetscObjectReference((PetscObject)mglevels[i]->A)); 556 } 557 } 558 559 PetscCall(KSPGetPC(mglevels[levels - 1]->smoothd, &tpc)); 560 PetscCall(PCPreSolveChangeRHS(tpc, &changed)); 561 PetscCall(KSPGetPC(mglevels[levels - 1]->smoothu, &tpc)); 562 PetscCall(PCPreSolveChangeRHS(tpc, &changeu)); 563 if (!changeu && !changed) { 564 if (matapp) { 565 PetscCall(MatDestroy(&mglevels[levels - 1]->B)); 566 mglevels[levels - 1]->B = B; 567 } else { 568 PetscCall(VecDestroy(&mglevels[levels - 1]->b)); 569 mglevels[levels - 1]->b = b; 570 } 571 } else { /* if the smoother changes the rhs during PreSolve, we cannot use the input vector */ 572 if (matapp) { 573 if (mglevels[levels - 1]->B) { 574 PetscInt N1, N2; 575 PetscBool flg; 576 577 PetscCall(MatGetSize(mglevels[levels - 1]->B, NULL, &N1)); 578 PetscCall(MatGetSize(B, NULL, &N2)); 579 PetscCall(PetscObjectTypeCompare((PetscObject)mglevels[levels - 1]->B, ((PetscObject)B)->type_name, &flg)); 580 if (N1 != N2 || !flg) PetscCall(MatDestroy(&mglevels[levels - 1]->B)); 581 } 582 if (!mglevels[levels - 1]->B) { 583 PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &mglevels[levels - 1]->B)); 584 } else { 585 PetscCall(MatCopy(B, mglevels[levels - 1]->B, SAME_NONZERO_PATTERN)); 586 } 587 } else { 588 if (!mglevels[levels - 1]->b) { 589 Vec *vec; 590 591 PetscCall(KSPCreateVecs(mglevels[levels - 1]->smoothd, 1, &vec, 0, NULL)); 592 mglevels[levels - 1]->b = *vec; 593 PetscCall(PetscFree(vec)); 594 } 595 PetscCall(VecCopy(b, mglevels[levels - 1]->b)); 596 } 597 } 598 if (matapp) { 599 mglevels[levels - 1]->X = X; 600 } else { 601 mglevels[levels - 1]->x = x; 602 } 603 604 /* If coarser Xs are present, it means we have already block applied the PC at least once 605 Reset operators if sizes/type do no match */ 606 if (matapp && levels > 1 && mglevels[levels - 2]->X) { 607 PetscInt Xc, Bc; 608 PetscBool flg; 609 610 PetscCall(MatGetSize(mglevels[levels - 2]->X, NULL, &Xc)); 611 PetscCall(MatGetSize(mglevels[levels - 1]->B, NULL, &Bc)); 612 PetscCall(PetscObjectTypeCompare((PetscObject)mglevels[levels - 2]->X, ((PetscObject)mglevels[levels - 1]->X)->type_name, &flg)); 613 if (Xc != Bc || !flg) { 614 PetscCall(MatDestroy(&mglevels[levels - 1]->R)); 615 for (i = 0; i < levels - 1; i++) { 616 PetscCall(MatDestroy(&mglevels[i]->R)); 617 PetscCall(MatDestroy(&mglevels[i]->B)); 618 PetscCall(MatDestroy(&mglevels[i]->X)); 619 } 620 } 621 } 622 623 if (mg->am == PC_MG_MULTIPLICATIVE) { 624 if (matapp) PetscCall(MatZeroEntries(X)); 625 else PetscCall(VecZeroEntries(x)); 626 for (i = 0; i < mg->cyclesperpcapply; i++) PetscCall(PCMGMCycle_Private(pc, mglevels + levels - 1, transpose, matapp, NULL)); 627 } else if (mg->am == PC_MG_ADDITIVE) { 628 PetscCall(PCMGACycle_Private(pc, mglevels, transpose, matapp)); 629 } else if (mg->am == PC_MG_KASKADE) { 630 PetscCall(PCMGKCycle_Private(pc, mglevels, transpose, matapp)); 631 } else { 632 PetscCall(PCMGFCycle_Private(pc, mglevels, transpose, matapp)); 633 } 634 if (mg->stageApply) PetscCall(PetscLogStagePop()); 635 if (!changeu && !changed) { 636 if (matapp) { 637 mglevels[levels - 1]->B = NULL; 638 } else { 639 mglevels[levels - 1]->b = NULL; 640 } 641 } 642 PetscFunctionReturn(PETSC_SUCCESS); 643 } 644 645 static PetscErrorCode PCApply_MG(PC pc, Vec b, Vec x) 646 { 647 PetscFunctionBegin; 648 PetscCall(PCApply_MG_Internal(pc, b, x, NULL, NULL, PETSC_FALSE)); 649 PetscFunctionReturn(PETSC_SUCCESS); 650 } 651 652 static PetscErrorCode PCApplyTranspose_MG(PC pc, Vec b, Vec x) 653 { 654 PetscFunctionBegin; 655 PetscCall(PCApply_MG_Internal(pc, b, x, NULL, NULL, PETSC_TRUE)); 656 PetscFunctionReturn(PETSC_SUCCESS); 657 } 658 659 static PetscErrorCode PCMatApply_MG(PC pc, Mat b, Mat x) 660 { 661 PetscFunctionBegin; 662 PetscCall(PCApply_MG_Internal(pc, NULL, NULL, b, x, PETSC_FALSE)); 663 PetscFunctionReturn(PETSC_SUCCESS); 664 } 665 666 PetscErrorCode PCSetFromOptions_MG(PC pc, PetscOptionItems *PetscOptionsObject) 667 { 668 PetscInt levels, cycles; 669 PetscBool flg, flg2; 670 PC_MG *mg = (PC_MG *)pc->data; 671 PC_MG_Levels **mglevels; 672 PCMGType mgtype; 673 PCMGCycleType mgctype; 674 PCMGGalerkinType gtype; 675 PCMGCoarseSpaceType coarseSpaceType; 676 677 PetscFunctionBegin; 678 levels = PetscMax(mg->nlevels, 1); 679 PetscOptionsHeadBegin(PetscOptionsObject, "Multigrid options"); 680 PetscCall(PetscOptionsInt("-pc_mg_levels", "Number of Levels", "PCMGSetLevels", levels, &levels, &flg)); 681 if (!flg && !mg->levels && pc->dm) { 682 PetscCall(DMGetRefineLevel(pc->dm, &levels)); 683 levels++; 684 mg->usedmfornumberoflevels = PETSC_TRUE; 685 } 686 PetscCall(PCMGSetLevels(pc, levels, NULL)); 687 mglevels = mg->levels; 688 689 mgctype = (PCMGCycleType)mglevels[0]->cycles; 690 PetscCall(PetscOptionsEnum("-pc_mg_cycle_type", "V cycle or for W-cycle", "PCMGSetCycleType", PCMGCycleTypes, (PetscEnum)mgctype, (PetscEnum *)&mgctype, &flg)); 691 if (flg) PetscCall(PCMGSetCycleType(pc, mgctype)); 692 coarseSpaceType = mg->coarseSpaceType; 693 PetscCall(PetscOptionsEnum("-pc_mg_adapt_interp_coarse_space", "Type of adaptive coarse space: none, polynomial, harmonic, eigenvector, generalized_eigenvector, gdsw", "PCMGSetAdaptCoarseSpaceType", PCMGCoarseSpaceTypes, (PetscEnum)coarseSpaceType, (PetscEnum *)&coarseSpaceType, &flg)); 694 if (flg) PetscCall(PCMGSetAdaptCoarseSpaceType(pc, coarseSpaceType)); 695 PetscCall(PetscOptionsInt("-pc_mg_adapt_interp_n", "Size of the coarse space for adaptive interpolation", "PCMGSetCoarseSpace", mg->Nc, &mg->Nc, &flg)); 696 PetscCall(PetscOptionsBool("-pc_mg_mesp_monitor", "Monitor the multilevel eigensolver", "PCMGSetAdaptInterpolation", PETSC_FALSE, &mg->mespMonitor, &flg)); 697 flg2 = PETSC_FALSE; 698 PetscCall(PetscOptionsBool("-pc_mg_adapt_cr", "Monitor coarse space quality using Compatible Relaxation (CR)", "PCMGSetAdaptCR", PETSC_FALSE, &flg2, &flg)); 699 if (flg) PetscCall(PCMGSetAdaptCR(pc, flg2)); 700 flg = PETSC_FALSE; 701 PetscCall(PetscOptionsBool("-pc_mg_distinct_smoothup", "Create separate smoothup KSP and append the prefix _up", "PCMGSetDistinctSmoothUp", PETSC_FALSE, &flg, NULL)); 702 if (flg) PetscCall(PCMGSetDistinctSmoothUp(pc)); 703 PetscCall(PetscOptionsEnum("-pc_mg_galerkin", "Use Galerkin process to compute coarser operators", "PCMGSetGalerkin", PCMGGalerkinTypes, (PetscEnum)mg->galerkin, (PetscEnum *)>ype, &flg)); 704 if (flg) PetscCall(PCMGSetGalerkin(pc, gtype)); 705 mgtype = mg->am; 706 PetscCall(PetscOptionsEnum("-pc_mg_type", "Multigrid type", "PCMGSetType", PCMGTypes, (PetscEnum)mgtype, (PetscEnum *)&mgtype, &flg)); 707 if (flg) PetscCall(PCMGSetType(pc, mgtype)); 708 if (mg->am == PC_MG_MULTIPLICATIVE) { 709 PetscCall(PetscOptionsInt("-pc_mg_multiplicative_cycles", "Number of cycles for each preconditioner step", "PCMGMultiplicativeSetCycles", mg->cyclesperpcapply, &cycles, &flg)); 710 if (flg) PetscCall(PCMGMultiplicativeSetCycles(pc, cycles)); 711 } 712 flg = PETSC_FALSE; 713 PetscCall(PetscOptionsBool("-pc_mg_log", "Log times for each multigrid level", "None", flg, &flg, NULL)); 714 if (flg) { 715 PetscInt i; 716 char eventname[128]; 717 718 levels = mglevels[0]->levels; 719 for (i = 0; i < levels; i++) { 720 PetscCall(PetscSNPrintf(eventname, PETSC_STATIC_ARRAY_LENGTH(eventname), "MGSetup Level %" PetscInt_FMT, i)); 721 PetscCall(PetscLogEventRegister(eventname, ((PetscObject)pc)->classid, &mglevels[i]->eventsmoothsetup)); 722 PetscCall(PetscSNPrintf(eventname, PETSC_STATIC_ARRAY_LENGTH(eventname), "MGSmooth Level %" PetscInt_FMT, i)); 723 PetscCall(PetscLogEventRegister(eventname, ((PetscObject)pc)->classid, &mglevels[i]->eventsmoothsolve)); 724 if (i) { 725 PetscCall(PetscSNPrintf(eventname, PETSC_STATIC_ARRAY_LENGTH(eventname), "MGResid Level %" PetscInt_FMT, i)); 726 PetscCall(PetscLogEventRegister(eventname, ((PetscObject)pc)->classid, &mglevels[i]->eventresidual)); 727 PetscCall(PetscSNPrintf(eventname, PETSC_STATIC_ARRAY_LENGTH(eventname), "MGInterp Level %" PetscInt_FMT, i)); 728 PetscCall(PetscLogEventRegister(eventname, ((PetscObject)pc)->classid, &mglevels[i]->eventinterprestrict)); 729 } 730 } 731 732 if (PetscDefined(USE_LOG)) { 733 const char sname[] = "MG Apply"; 734 735 PetscCall(PetscLogStageGetId(sname, &mg->stageApply)); 736 if (mg->stageApply < 0) PetscCall(PetscLogStageRegister(sname, &mg->stageApply)); 737 } 738 } 739 PetscOptionsHeadEnd(); 740 /* Check option consistency */ 741 PetscCall(PCMGGetGalerkin(pc, >ype)); 742 PetscCall(PCMGGetAdaptInterpolation(pc, &flg)); 743 PetscCheck(!flg || !(gtype >= PC_MG_GALERKIN_NONE), PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Must use Galerkin coarse operators when adapting the interpolator"); 744 PetscFunctionReturn(PETSC_SUCCESS); 745 } 746 747 const char *const PCMGTypes[] = {"MULTIPLICATIVE", "ADDITIVE", "FULL", "KASKADE", "PCMGType", "PC_MG", NULL}; 748 const char *const PCMGCycleTypes[] = {"invalid", "v", "w", "PCMGCycleType", "PC_MG_CYCLE", NULL}; 749 const char *const PCMGGalerkinTypes[] = {"both", "pmat", "mat", "none", "external", "PCMGGalerkinType", "PC_MG_GALERKIN", NULL}; 750 const char *const PCMGCoarseSpaceTypes[] = {"none", "polynomial", "harmonic", "eigenvector", "generalized_eigenvector", "gdsw", "PCMGCoarseSpaceType", "PCMG_ADAPT_NONE", NULL}; 751 752 #include <petscdraw.h> 753 PetscErrorCode PCView_MG(PC pc, PetscViewer viewer) 754 { 755 PC_MG *mg = (PC_MG *)pc->data; 756 PC_MG_Levels **mglevels = mg->levels; 757 PetscInt levels = mglevels ? mglevels[0]->levels : 0, i; 758 PetscBool iascii, isbinary, isdraw; 759 760 PetscFunctionBegin; 761 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 762 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary)); 763 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw)); 764 if (iascii) { 765 const char *cyclename = levels ? (mglevels[0]->cycles == PC_MG_CYCLE_V ? "v" : "w") : "unknown"; 766 PetscCall(PetscViewerASCIIPrintf(viewer, " type is %s, levels=%" PetscInt_FMT " cycles=%s\n", PCMGTypes[mg->am], levels, cyclename)); 767 if (mg->am == PC_MG_MULTIPLICATIVE) PetscCall(PetscViewerASCIIPrintf(viewer, " Cycles per PCApply=%" PetscInt_FMT "\n", mg->cyclesperpcapply)); 768 if (mg->galerkin == PC_MG_GALERKIN_BOTH) { 769 PetscCall(PetscViewerASCIIPrintf(viewer, " Using Galerkin computed coarse grid matrices\n")); 770 } else if (mg->galerkin == PC_MG_GALERKIN_PMAT) { 771 PetscCall(PetscViewerASCIIPrintf(viewer, " Using Galerkin computed coarse grid matrices for pmat\n")); 772 } else if (mg->galerkin == PC_MG_GALERKIN_MAT) { 773 PetscCall(PetscViewerASCIIPrintf(viewer, " Using Galerkin computed coarse grid matrices for mat\n")); 774 } else if (mg->galerkin == PC_MG_GALERKIN_EXTERNAL) { 775 PetscCall(PetscViewerASCIIPrintf(viewer, " Using externally compute Galerkin coarse grid matrices\n")); 776 } else { 777 PetscCall(PetscViewerASCIIPrintf(viewer, " Not using Galerkin computed coarse grid matrices\n")); 778 } 779 if (mg->view) PetscCall((*mg->view)(pc, viewer)); 780 for (i = 0; i < levels; i++) { 781 if (i) { 782 PetscCall(PetscViewerASCIIPrintf(viewer, "Down solver (pre-smoother) on level %" PetscInt_FMT " -------------------------------\n", i)); 783 } else { 784 PetscCall(PetscViewerASCIIPrintf(viewer, "Coarse grid solver -- level %" PetscInt_FMT " -------------------------------\n", i)); 785 } 786 PetscCall(PetscViewerASCIIPushTab(viewer)); 787 PetscCall(KSPView(mglevels[i]->smoothd, viewer)); 788 PetscCall(PetscViewerASCIIPopTab(viewer)); 789 if (i && mglevels[i]->smoothd == mglevels[i]->smoothu) { 790 PetscCall(PetscViewerASCIIPrintf(viewer, "Up solver (post-smoother) same as down solver (pre-smoother)\n")); 791 } else if (i) { 792 PetscCall(PetscViewerASCIIPrintf(viewer, "Up solver (post-smoother) on level %" PetscInt_FMT " -------------------------------\n", i)); 793 PetscCall(PetscViewerASCIIPushTab(viewer)); 794 PetscCall(KSPView(mglevels[i]->smoothu, viewer)); 795 PetscCall(PetscViewerASCIIPopTab(viewer)); 796 } 797 if (i && mglevels[i]->cr) { 798 PetscCall(PetscViewerASCIIPrintf(viewer, "CR solver on level %" PetscInt_FMT " -------------------------------\n", i)); 799 PetscCall(PetscViewerASCIIPushTab(viewer)); 800 PetscCall(KSPView(mglevels[i]->cr, viewer)); 801 PetscCall(PetscViewerASCIIPopTab(viewer)); 802 } 803 } 804 } else if (isbinary) { 805 for (i = levels - 1; i >= 0; i--) { 806 PetscCall(KSPView(mglevels[i]->smoothd, viewer)); 807 if (i && mglevels[i]->smoothd != mglevels[i]->smoothu) PetscCall(KSPView(mglevels[i]->smoothu, viewer)); 808 } 809 } else if (isdraw) { 810 PetscDraw draw; 811 PetscReal x, w, y, bottom, th; 812 PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw)); 813 PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y)); 814 PetscCall(PetscDrawStringGetSize(draw, NULL, &th)); 815 bottom = y - th; 816 for (i = levels - 1; i >= 0; i--) { 817 if (!mglevels[i]->smoothu || (mglevels[i]->smoothu == mglevels[i]->smoothd)) { 818 PetscCall(PetscDrawPushCurrentPoint(draw, x, bottom)); 819 PetscCall(KSPView(mglevels[i]->smoothd, viewer)); 820 PetscCall(PetscDrawPopCurrentPoint(draw)); 821 } else { 822 w = 0.5 * PetscMin(1.0 - x, x); 823 PetscCall(PetscDrawPushCurrentPoint(draw, x + w, bottom)); 824 PetscCall(KSPView(mglevels[i]->smoothd, viewer)); 825 PetscCall(PetscDrawPopCurrentPoint(draw)); 826 PetscCall(PetscDrawPushCurrentPoint(draw, x - w, bottom)); 827 PetscCall(KSPView(mglevels[i]->smoothu, viewer)); 828 PetscCall(PetscDrawPopCurrentPoint(draw)); 829 } 830 PetscCall(PetscDrawGetBoundingBox(draw, NULL, &bottom, NULL, NULL)); 831 bottom -= th; 832 } 833 } 834 PetscFunctionReturn(PETSC_SUCCESS); 835 } 836 837 #include <petsc/private/kspimpl.h> 838 839 /* 840 Calls setup for the KSP on each level 841 */ 842 PetscErrorCode PCSetUp_MG(PC pc) 843 { 844 PC_MG *mg = (PC_MG *)pc->data; 845 PC_MG_Levels **mglevels = mg->levels; 846 PetscInt i, n; 847 PC cpc; 848 PetscBool dump = PETSC_FALSE, opsset, use_amat, missinginterpolate = PETSC_FALSE; 849 Mat dA, dB; 850 Vec tvec; 851 DM *dms; 852 PetscViewer viewer = NULL; 853 PetscBool dAeqdB = PETSC_FALSE, needRestricts = PETSC_FALSE, doCR = PETSC_FALSE; 854 PetscBool adaptInterpolation = mg->adaptInterpolation; 855 856 PetscFunctionBegin; 857 PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set MG levels with PCMGSetLevels() before setting up"); 858 n = mglevels[0]->levels; 859 /* FIX: Move this to PCSetFromOptions_MG? */ 860 if (mg->usedmfornumberoflevels) { 861 PetscInt levels; 862 PetscCall(DMGetRefineLevel(pc->dm, &levels)); 863 levels++; 864 if (levels > n) { /* the problem is now being solved on a finer grid */ 865 PetscCall(PCMGSetLevels(pc, levels, NULL)); 866 n = levels; 867 PetscCall(PCSetFromOptions(pc)); /* it is bad to call this here, but otherwise will never be called for the new hierarchy */ 868 mglevels = mg->levels; 869 } 870 } 871 PetscCall(KSPGetPC(mglevels[0]->smoothd, &cpc)); 872 873 /* If user did not provide fine grid operators OR operator was not updated since last global KSPSetOperators() */ 874 /* so use those from global PC */ 875 /* Is this what we always want? What if user wants to keep old one? */ 876 PetscCall(KSPGetOperatorsSet(mglevels[n - 1]->smoothd, NULL, &opsset)); 877 if (opsset) { 878 Mat mmat; 879 PetscCall(KSPGetOperators(mglevels[n - 1]->smoothd, NULL, &mmat)); 880 if (mmat == pc->pmat) opsset = PETSC_FALSE; 881 } 882 883 /* Create CR solvers */ 884 PetscCall(PCMGGetAdaptCR(pc, &doCR)); 885 if (doCR) { 886 const char *prefix; 887 888 PetscCall(PCGetOptionsPrefix(pc, &prefix)); 889 for (i = 1; i < n; ++i) { 890 PC ipc, cr; 891 char crprefix[128]; 892 893 PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &mglevels[i]->cr)); 894 PetscCall(KSPSetNestLevel(mglevels[i]->cr, pc->kspnestlevel)); 895 PetscCall(KSPSetErrorIfNotConverged(mglevels[i]->cr, PETSC_FALSE)); 896 PetscCall(PetscObjectIncrementTabLevel((PetscObject)mglevels[i]->cr, (PetscObject)pc, n - i)); 897 PetscCall(KSPSetOptionsPrefix(mglevels[i]->cr, prefix)); 898 PetscCall(PetscObjectComposedDataSetInt((PetscObject)mglevels[i]->cr, PetscMGLevelId, mglevels[i]->level)); 899 PetscCall(KSPSetType(mglevels[i]->cr, KSPCHEBYSHEV)); 900 PetscCall(KSPSetConvergenceTest(mglevels[i]->cr, KSPConvergedSkip, NULL, NULL)); 901 PetscCall(KSPSetNormType(mglevels[i]->cr, KSP_NORM_PRECONDITIONED)); 902 PetscCall(KSPGetPC(mglevels[i]->cr, &ipc)); 903 904 PetscCall(PCSetType(ipc, PCCOMPOSITE)); 905 PetscCall(PCCompositeSetType(ipc, PC_COMPOSITE_MULTIPLICATIVE)); 906 PetscCall(PCCompositeAddPCType(ipc, PCSOR)); 907 PetscCall(CreateCR_Private(pc, i, &cr)); 908 PetscCall(PCCompositeAddPC(ipc, cr)); 909 PetscCall(PCDestroy(&cr)); 910 911 PetscCall(KSPSetTolerances(mglevels[i]->cr, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT, mg->default_smoothd)); 912 PetscCall(KSPSetInitialGuessNonzero(mglevels[i]->cr, PETSC_TRUE)); 913 PetscCall(PetscSNPrintf(crprefix, 128, "mg_levels_%" PetscInt_FMT "_cr_", i)); 914 PetscCall(KSPAppendOptionsPrefix(mglevels[i]->cr, crprefix)); 915 } 916 } 917 918 if (!opsset) { 919 PetscCall(PCGetUseAmat(pc, &use_amat)); 920 if (use_amat) { 921 PetscCall(PetscInfo(pc, "Using outer operators to define finest grid operator \n because PCMGGetSmoother(pc,nlevels-1,&ksp);KSPSetOperators(ksp,...); was not called.\n")); 922 PetscCall(KSPSetOperators(mglevels[n - 1]->smoothd, pc->mat, pc->pmat)); 923 } else { 924 PetscCall(PetscInfo(pc, "Using matrix (pmat) operators to define finest grid operator \n because PCMGGetSmoother(pc,nlevels-1,&ksp);KSPSetOperators(ksp,...); was not called.\n")); 925 PetscCall(KSPSetOperators(mglevels[n - 1]->smoothd, pc->pmat, pc->pmat)); 926 } 927 } 928 929 for (i = n - 1; i > 0; i--) { 930 if (!(mglevels[i]->interpolate || mglevels[i]->restrct)) { 931 missinginterpolate = PETSC_TRUE; 932 break; 933 } 934 } 935 936 PetscCall(KSPGetOperators(mglevels[n - 1]->smoothd, &dA, &dB)); 937 if (dA == dB) dAeqdB = PETSC_TRUE; 938 if (mg->galerkin == PC_MG_GALERKIN_NONE || ((mg->galerkin == PC_MG_GALERKIN_PMAT || mg->galerkin == PC_MG_GALERKIN_MAT) && !dAeqdB)) { 939 needRestricts = PETSC_TRUE; /* user must compute either mat, pmat, or both so must restrict x to coarser levels */ 940 } 941 942 if (pc->dm && !pc->setupcalled) { 943 /* finest smoother also gets DM but it is not active, independent of whether galerkin==PC_MG_GALERKIN_EXTERNAL */ 944 PetscCall(KSPSetDM(mglevels[n - 1]->smoothd, pc->dm)); 945 PetscCall(KSPSetDMActive(mglevels[n - 1]->smoothd, PETSC_FALSE)); 946 if (mglevels[n - 1]->smoothd != mglevels[n - 1]->smoothu) { 947 PetscCall(KSPSetDM(mglevels[n - 1]->smoothu, pc->dm)); 948 PetscCall(KSPSetDMActive(mglevels[n - 1]->smoothu, PETSC_FALSE)); 949 } 950 if (mglevels[n - 1]->cr) { 951 PetscCall(KSPSetDM(mglevels[n - 1]->cr, pc->dm)); 952 PetscCall(KSPSetDMActive(mglevels[n - 1]->cr, PETSC_FALSE)); 953 } 954 } 955 956 /* 957 Skipping if user has provided all interpolation/restriction needed (since DM might not be able to produce them (when coming from SNES/TS) 958 Skipping for externally managed hierarchy (such as ML and GAMG). Cleaner logic here would be great. Wrap ML/GAMG as DMs? 959 */ 960 if (missinginterpolate && mg->galerkin != PC_MG_GALERKIN_EXTERNAL && !pc->setupcalled) { 961 /* first see if we can compute a coarse space */ 962 if (mg->coarseSpaceType == PCMG_ADAPT_GDSW) { 963 for (i = n - 2; i > -1; i--) { 964 if (!mglevels[i + 1]->restrct && !mglevels[i + 1]->interpolate) { 965 PetscCall(PCMGComputeCoarseSpace_Internal(pc, i + 1, mg->coarseSpaceType, mg->Nc, NULL, &mglevels[i + 1]->coarseSpace)); 966 PetscCall(PCMGSetInterpolation(pc, i + 1, mglevels[i + 1]->coarseSpace)); 967 } 968 } 969 } else { /* construct the interpolation from the DMs */ 970 Mat p; 971 Vec rscale; 972 PetscCall(PetscMalloc1(n, &dms)); 973 dms[n - 1] = pc->dm; 974 /* Separately create them so we do not get DMKSP interference between levels */ 975 for (i = n - 2; i > -1; i--) PetscCall(DMCoarsen(dms[i + 1], MPI_COMM_NULL, &dms[i])); 976 for (i = n - 2; i > -1; i--) { 977 DMKSP kdm; 978 PetscBool dmhasrestrict, dmhasinject; 979 PetscCall(KSPSetDM(mglevels[i]->smoothd, dms[i])); 980 if (!needRestricts) PetscCall(KSPSetDMActive(mglevels[i]->smoothd, PETSC_FALSE)); 981 if (mglevels[i]->smoothd != mglevels[i]->smoothu) { 982 PetscCall(KSPSetDM(mglevels[i]->smoothu, dms[i])); 983 if (!needRestricts) PetscCall(KSPSetDMActive(mglevels[i]->smoothu, PETSC_FALSE)); 984 } 985 if (mglevels[i]->cr) { 986 PetscCall(KSPSetDM(mglevels[i]->cr, dms[i])); 987 if (!needRestricts) PetscCall(KSPSetDMActive(mglevels[i]->cr, PETSC_FALSE)); 988 } 989 PetscCall(DMGetDMKSPWrite(dms[i], &kdm)); 990 /* Ugly hack so that the next KSPSetUp() will use the RHS that we set. A better fix is to change dmActive to take 991 * a bitwise OR of computing the matrix, RHS, and initial iterate. */ 992 kdm->ops->computerhs = NULL; 993 kdm->rhsctx = NULL; 994 if (!mglevels[i + 1]->interpolate) { 995 PetscCall(DMCreateInterpolation(dms[i], dms[i + 1], &p, &rscale)); 996 PetscCall(PCMGSetInterpolation(pc, i + 1, p)); 997 if (rscale) PetscCall(PCMGSetRScale(pc, i + 1, rscale)); 998 PetscCall(VecDestroy(&rscale)); 999 PetscCall(MatDestroy(&p)); 1000 } 1001 PetscCall(DMHasCreateRestriction(dms[i], &dmhasrestrict)); 1002 if (dmhasrestrict && !mglevels[i + 1]->restrct) { 1003 PetscCall(DMCreateRestriction(dms[i], dms[i + 1], &p)); 1004 PetscCall(PCMGSetRestriction(pc, i + 1, p)); 1005 PetscCall(MatDestroy(&p)); 1006 } 1007 PetscCall(DMHasCreateInjection(dms[i], &dmhasinject)); 1008 if (dmhasinject && !mglevels[i + 1]->inject) { 1009 PetscCall(DMCreateInjection(dms[i], dms[i + 1], &p)); 1010 PetscCall(PCMGSetInjection(pc, i + 1, p)); 1011 PetscCall(MatDestroy(&p)); 1012 } 1013 } 1014 1015 for (i = n - 2; i > -1; i--) PetscCall(DMDestroy(&dms[i])); 1016 PetscCall(PetscFree(dms)); 1017 } 1018 } 1019 1020 if (mg->galerkin < PC_MG_GALERKIN_NONE) { 1021 Mat A, B; 1022 PetscBool doA = PETSC_FALSE, doB = PETSC_FALSE; 1023 MatReuse reuse = MAT_INITIAL_MATRIX; 1024 1025 if (mg->galerkin == PC_MG_GALERKIN_PMAT || mg->galerkin == PC_MG_GALERKIN_BOTH) doB = PETSC_TRUE; 1026 if (mg->galerkin == PC_MG_GALERKIN_MAT || (mg->galerkin == PC_MG_GALERKIN_BOTH && dA != dB)) doA = PETSC_TRUE; 1027 if (pc->setupcalled) reuse = MAT_REUSE_MATRIX; 1028 for (i = n - 2; i > -1; i--) { 1029 PetscCheck(mglevels[i + 1]->restrct || mglevels[i + 1]->interpolate, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must provide interpolation or restriction for each MG level except level 0"); 1030 if (!mglevels[i + 1]->interpolate) PetscCall(PCMGSetInterpolation(pc, i + 1, mglevels[i + 1]->restrct)); 1031 if (!mglevels[i + 1]->restrct) PetscCall(PCMGSetRestriction(pc, i + 1, mglevels[i + 1]->interpolate)); 1032 if (reuse == MAT_REUSE_MATRIX) PetscCall(KSPGetOperators(mglevels[i]->smoothd, &A, &B)); 1033 if (doA) PetscCall(MatGalerkin(mglevels[i + 1]->restrct, dA, mglevels[i + 1]->interpolate, reuse, 1.0, &A)); 1034 if (doB) PetscCall(MatGalerkin(mglevels[i + 1]->restrct, dB, mglevels[i + 1]->interpolate, reuse, 1.0, &B)); 1035 /* the management of the PetscObjectReference() and PetscObjecDereference() below is rather delicate */ 1036 if (!doA && dAeqdB) { 1037 if (reuse == MAT_INITIAL_MATRIX) PetscCall(PetscObjectReference((PetscObject)B)); 1038 A = B; 1039 } else if (!doA && reuse == MAT_INITIAL_MATRIX) { 1040 PetscCall(KSPGetOperators(mglevels[i]->smoothd, &A, NULL)); 1041 PetscCall(PetscObjectReference((PetscObject)A)); 1042 } 1043 if (!doB && dAeqdB) { 1044 if (reuse == MAT_INITIAL_MATRIX) PetscCall(PetscObjectReference((PetscObject)A)); 1045 B = A; 1046 } else if (!doB && reuse == MAT_INITIAL_MATRIX) { 1047 PetscCall(KSPGetOperators(mglevels[i]->smoothd, NULL, &B)); 1048 PetscCall(PetscObjectReference((PetscObject)B)); 1049 } 1050 if (reuse == MAT_INITIAL_MATRIX) { 1051 PetscCall(KSPSetOperators(mglevels[i]->smoothd, A, B)); 1052 PetscCall(PetscObjectDereference((PetscObject)A)); 1053 PetscCall(PetscObjectDereference((PetscObject)B)); 1054 } 1055 dA = A; 1056 dB = B; 1057 } 1058 } 1059 1060 /* Adapt interpolation matrices */ 1061 if (adaptInterpolation) { 1062 for (i = 0; i < n; ++i) { 1063 if (!mglevels[i]->coarseSpace) PetscCall(PCMGComputeCoarseSpace_Internal(pc, i, mg->coarseSpaceType, mg->Nc, !i ? NULL : mglevels[i - 1]->coarseSpace, &mglevels[i]->coarseSpace)); 1064 if (i) PetscCall(PCMGAdaptInterpolator_Internal(pc, i, mglevels[i - 1]->smoothu, mglevels[i]->smoothu, mglevels[i - 1]->coarseSpace, mglevels[i]->coarseSpace)); 1065 } 1066 for (i = n - 2; i > -1; --i) PetscCall(PCMGRecomputeLevelOperators_Internal(pc, i)); 1067 } 1068 1069 if (needRestricts && pc->dm) { 1070 for (i = n - 2; i >= 0; i--) { 1071 DM dmfine, dmcoarse; 1072 Mat Restrict, Inject; 1073 Vec rscale; 1074 PetscCall(KSPGetDM(mglevels[i + 1]->smoothd, &dmfine)); 1075 PetscCall(KSPGetDM(mglevels[i]->smoothd, &dmcoarse)); 1076 PetscCall(PCMGGetRestriction(pc, i + 1, &Restrict)); 1077 PetscCall(PCMGGetRScale(pc, i + 1, &rscale)); 1078 PetscCall(PCMGGetInjection(pc, i + 1, &Inject)); 1079 PetscCall(DMRestrict(dmfine, Restrict, rscale, Inject, dmcoarse)); 1080 } 1081 } 1082 1083 if (!pc->setupcalled) { 1084 for (i = 0; i < n; i++) PetscCall(KSPSetFromOptions(mglevels[i]->smoothd)); 1085 for (i = 1; i < n; i++) { 1086 if (mglevels[i]->smoothu && (mglevels[i]->smoothu != mglevels[i]->smoothd)) PetscCall(KSPSetFromOptions(mglevels[i]->smoothu)); 1087 if (mglevels[i]->cr) PetscCall(KSPSetFromOptions(mglevels[i]->cr)); 1088 } 1089 /* insure that if either interpolation or restriction is set the other one is set */ 1090 for (i = 1; i < n; i++) { 1091 PetscCall(PCMGGetInterpolation(pc, i, NULL)); 1092 PetscCall(PCMGGetRestriction(pc, i, NULL)); 1093 } 1094 for (i = 0; i < n - 1; i++) { 1095 if (!mglevels[i]->b) { 1096 Vec *vec; 1097 PetscCall(KSPCreateVecs(mglevels[i]->smoothd, 1, &vec, 0, NULL)); 1098 PetscCall(PCMGSetRhs(pc, i, *vec)); 1099 PetscCall(VecDestroy(vec)); 1100 PetscCall(PetscFree(vec)); 1101 } 1102 if (!mglevels[i]->r && i) { 1103 PetscCall(VecDuplicate(mglevels[i]->b, &tvec)); 1104 PetscCall(PCMGSetR(pc, i, tvec)); 1105 PetscCall(VecDestroy(&tvec)); 1106 } 1107 if (!mglevels[i]->x) { 1108 PetscCall(VecDuplicate(mglevels[i]->b, &tvec)); 1109 PetscCall(PCMGSetX(pc, i, tvec)); 1110 PetscCall(VecDestroy(&tvec)); 1111 } 1112 if (doCR) { 1113 PetscCall(VecDuplicate(mglevels[i]->b, &mglevels[i]->crx)); 1114 PetscCall(VecDuplicate(mglevels[i]->b, &mglevels[i]->crb)); 1115 PetscCall(VecZeroEntries(mglevels[i]->crb)); 1116 } 1117 } 1118 if (n != 1 && !mglevels[n - 1]->r) { 1119 /* PCMGSetR() on the finest level if user did not supply it */ 1120 Vec *vec; 1121 PetscCall(KSPCreateVecs(mglevels[n - 1]->smoothd, 1, &vec, 0, NULL)); 1122 PetscCall(PCMGSetR(pc, n - 1, *vec)); 1123 PetscCall(VecDestroy(vec)); 1124 PetscCall(PetscFree(vec)); 1125 } 1126 if (doCR) { 1127 PetscCall(VecDuplicate(mglevels[n - 1]->r, &mglevels[n - 1]->crx)); 1128 PetscCall(VecDuplicate(mglevels[n - 1]->r, &mglevels[n - 1]->crb)); 1129 PetscCall(VecZeroEntries(mglevels[n - 1]->crb)); 1130 } 1131 } 1132 1133 if (pc->dm) { 1134 /* need to tell all the coarser levels to rebuild the matrix using the DM for that level */ 1135 for (i = 0; i < n - 1; i++) { 1136 if (mglevels[i]->smoothd->setupstage != KSP_SETUP_NEW) mglevels[i]->smoothd->setupstage = KSP_SETUP_NEWMATRIX; 1137 } 1138 } 1139 // We got here (PCSetUp_MG) because the matrix has changed, which means the smoother needs to be set up again (e.g., 1140 // new diagonal for Jacobi). Setting it here allows it to be logged under PCSetUp rather than deep inside a PCApply. 1141 if (mglevels[n - 1]->smoothd->setupstage != KSP_SETUP_NEW) mglevels[n - 1]->smoothd->setupstage = KSP_SETUP_NEWMATRIX; 1142 1143 for (i = 1; i < n; i++) { 1144 if (mglevels[i]->smoothu == mglevels[i]->smoothd || mg->am == PC_MG_FULL || mg->am == PC_MG_KASKADE || mg->cyclesperpcapply > 1) { 1145 /* if doing only down then initial guess is zero */ 1146 PetscCall(KSPSetInitialGuessNonzero(mglevels[i]->smoothd, PETSC_TRUE)); 1147 } 1148 if (mglevels[i]->cr) PetscCall(KSPSetInitialGuessNonzero(mglevels[i]->cr, PETSC_TRUE)); 1149 if (mglevels[i]->eventsmoothsetup) PetscCall(PetscLogEventBegin(mglevels[i]->eventsmoothsetup, 0, 0, 0, 0)); 1150 PetscCall(KSPSetUp(mglevels[i]->smoothd)); 1151 if (mglevels[i]->smoothd->reason) pc->failedreason = PC_SUBPC_ERROR; 1152 if (mglevels[i]->eventsmoothsetup) PetscCall(PetscLogEventEnd(mglevels[i]->eventsmoothsetup, 0, 0, 0, 0)); 1153 if (!mglevels[i]->residual) { 1154 Mat mat; 1155 PetscCall(KSPGetOperators(mglevels[i]->smoothd, &mat, NULL)); 1156 PetscCall(PCMGSetResidual(pc, i, PCMGResidualDefault, mat)); 1157 } 1158 if (!mglevels[i]->residualtranspose) { 1159 Mat mat; 1160 PetscCall(KSPGetOperators(mglevels[i]->smoothd, &mat, NULL)); 1161 PetscCall(PCMGSetResidualTranspose(pc, i, PCMGResidualTransposeDefault, mat)); 1162 } 1163 } 1164 for (i = 1; i < n; i++) { 1165 if (mglevels[i]->smoothu && mglevels[i]->smoothu != mglevels[i]->smoothd) { 1166 Mat downmat, downpmat; 1167 1168 /* check if operators have been set for up, if not use down operators to set them */ 1169 PetscCall(KSPGetOperatorsSet(mglevels[i]->smoothu, &opsset, NULL)); 1170 if (!opsset) { 1171 PetscCall(KSPGetOperators(mglevels[i]->smoothd, &downmat, &downpmat)); 1172 PetscCall(KSPSetOperators(mglevels[i]->smoothu, downmat, downpmat)); 1173 } 1174 1175 PetscCall(KSPSetInitialGuessNonzero(mglevels[i]->smoothu, PETSC_TRUE)); 1176 if (mglevels[i]->eventsmoothsetup) PetscCall(PetscLogEventBegin(mglevels[i]->eventsmoothsetup, 0, 0, 0, 0)); 1177 PetscCall(KSPSetUp(mglevels[i]->smoothu)); 1178 if (mglevels[i]->smoothu->reason) pc->failedreason = PC_SUBPC_ERROR; 1179 if (mglevels[i]->eventsmoothsetup) PetscCall(PetscLogEventEnd(mglevels[i]->eventsmoothsetup, 0, 0, 0, 0)); 1180 } 1181 if (mglevels[i]->cr) { 1182 Mat downmat, downpmat; 1183 1184 /* check if operators have been set for up, if not use down operators to set them */ 1185 PetscCall(KSPGetOperatorsSet(mglevels[i]->cr, &opsset, NULL)); 1186 if (!opsset) { 1187 PetscCall(KSPGetOperators(mglevels[i]->smoothd, &downmat, &downpmat)); 1188 PetscCall(KSPSetOperators(mglevels[i]->cr, downmat, downpmat)); 1189 } 1190 1191 PetscCall(KSPSetInitialGuessNonzero(mglevels[i]->cr, PETSC_TRUE)); 1192 if (mglevels[i]->eventsmoothsetup) PetscCall(PetscLogEventBegin(mglevels[i]->eventsmoothsetup, 0, 0, 0, 0)); 1193 PetscCall(KSPSetUp(mglevels[i]->cr)); 1194 if (mglevels[i]->cr->reason) pc->failedreason = PC_SUBPC_ERROR; 1195 if (mglevels[i]->eventsmoothsetup) PetscCall(PetscLogEventEnd(mglevels[i]->eventsmoothsetup, 0, 0, 0, 0)); 1196 } 1197 } 1198 1199 if (mglevels[0]->eventsmoothsetup) PetscCall(PetscLogEventBegin(mglevels[0]->eventsmoothsetup, 0, 0, 0, 0)); 1200 PetscCall(KSPSetUp(mglevels[0]->smoothd)); 1201 if (mglevels[0]->smoothd->reason) pc->failedreason = PC_SUBPC_ERROR; 1202 if (mglevels[0]->eventsmoothsetup) PetscCall(PetscLogEventEnd(mglevels[0]->eventsmoothsetup, 0, 0, 0, 0)); 1203 1204 /* 1205 Dump the interpolation/restriction matrices plus the 1206 Jacobian/stiffness on each level. This allows MATLAB users to 1207 easily check if the Galerkin condition A_c = R A_f R^T is satisfied. 1208 1209 Only support one or the other at the same time. 1210 */ 1211 #if defined(PETSC_USE_SOCKET_VIEWER) 1212 PetscCall(PetscOptionsGetBool(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_mg_dump_matlab", &dump, NULL)); 1213 if (dump) viewer = PETSC_VIEWER_SOCKET_(PetscObjectComm((PetscObject)pc)); 1214 dump = PETSC_FALSE; 1215 #endif 1216 PetscCall(PetscOptionsGetBool(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_mg_dump_binary", &dump, NULL)); 1217 if (dump) viewer = PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)pc)); 1218 1219 if (viewer) { 1220 for (i = 1; i < n; i++) PetscCall(MatView(mglevels[i]->restrct, viewer)); 1221 for (i = 0; i < n; i++) { 1222 PetscCall(KSPGetPC(mglevels[i]->smoothd, &pc)); 1223 PetscCall(MatView(pc->mat, viewer)); 1224 } 1225 } 1226 PetscFunctionReturn(PETSC_SUCCESS); 1227 } 1228 1229 PetscErrorCode PCMGGetLevels_MG(PC pc, PetscInt *levels) 1230 { 1231 PC_MG *mg = (PC_MG *)pc->data; 1232 1233 PetscFunctionBegin; 1234 *levels = mg->nlevels; 1235 PetscFunctionReturn(PETSC_SUCCESS); 1236 } 1237 1238 /*@ 1239 PCMGGetLevels - Gets the number of levels to use with `PCMG`. 1240 1241 Not Collective 1242 1243 Input Parameter: 1244 . pc - the preconditioner context 1245 1246 Output Parameter: 1247 . levels - the number of levels 1248 1249 Level: advanced 1250 1251 .seealso: [](ch_ksp), `PCMG`, `PCMGSetLevels()` 1252 @*/ 1253 PetscErrorCode PCMGGetLevels(PC pc, PetscInt *levels) 1254 { 1255 PetscFunctionBegin; 1256 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1257 PetscAssertPointer(levels, 2); 1258 *levels = 0; 1259 PetscTryMethod(pc, "PCMGGetLevels_C", (PC, PetscInt *), (pc, levels)); 1260 PetscFunctionReturn(PETSC_SUCCESS); 1261 } 1262 1263 /*@ 1264 PCMGGetGridComplexity - compute operator and grid complexity of the `PCMG` hierarchy 1265 1266 Input Parameter: 1267 . pc - the preconditioner context 1268 1269 Output Parameters: 1270 + gc - grid complexity = sum_i(n_i) / n_0 1271 - oc - operator complexity = sum_i(nnz_i) / nnz_0 1272 1273 Level: advanced 1274 1275 Note: 1276 This is often call the operator complexity in multigrid literature 1277 1278 .seealso: [](ch_ksp), `PCMG`, `PCMGGetLevels()`, `PCMGSetLevels()` 1279 @*/ 1280 PetscErrorCode PCMGGetGridComplexity(PC pc, PetscReal *gc, PetscReal *oc) 1281 { 1282 PC_MG *mg = (PC_MG *)pc->data; 1283 PC_MG_Levels **mglevels = mg->levels; 1284 PetscInt lev, N; 1285 PetscLogDouble nnz0 = 0, sgc = 0, soc = 0, n0 = 0; 1286 MatInfo info; 1287 1288 PetscFunctionBegin; 1289 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1290 if (gc) PetscAssertPointer(gc, 2); 1291 if (oc) PetscAssertPointer(oc, 3); 1292 if (!pc->setupcalled) { 1293 if (gc) *gc = 0; 1294 if (oc) *oc = 0; 1295 PetscFunctionReturn(PETSC_SUCCESS); 1296 } 1297 PetscCheck(mg->nlevels > 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MG has no levels"); 1298 for (lev = 0; lev < mg->nlevels; lev++) { 1299 Mat dB; 1300 PetscCall(KSPGetOperators(mglevels[lev]->smoothd, NULL, &dB)); 1301 PetscCall(MatGetInfo(dB, MAT_GLOBAL_SUM, &info)); /* global reduction */ 1302 PetscCall(MatGetSize(dB, &N, NULL)); 1303 sgc += N; 1304 soc += info.nz_used; 1305 if (lev == mg->nlevels - 1) { 1306 nnz0 = info.nz_used; 1307 n0 = N; 1308 } 1309 } 1310 PetscCheck(n0 > 0 && gc, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number for grid points on finest level is not available"); 1311 *gc = (PetscReal)(sgc / n0); 1312 if (nnz0 > 0 && oc) *oc = (PetscReal)(soc / nnz0); 1313 PetscFunctionReturn(PETSC_SUCCESS); 1314 } 1315 1316 /*@ 1317 PCMGSetType - Determines the form of multigrid to use, either 1318 multiplicative, additive, full, or the Kaskade algorithm. 1319 1320 Logically Collective 1321 1322 Input Parameters: 1323 + pc - the preconditioner context 1324 - form - multigrid form, one of `PC_MG_MULTIPLICATIVE`, `PC_MG_ADDITIVE`, `PC_MG_FULL`, `PC_MG_KASKADE` 1325 1326 Options Database Key: 1327 . -pc_mg_type <form> - Sets <form>, one of multiplicative, additive, full, kaskade 1328 1329 Level: advanced 1330 1331 .seealso: [](ch_ksp), `PCMGType`, `PCMG`, `PCMGGetLevels()`, `PCMGSetLevels()`, `PCMGGetType()`, `PCMGCycleType` 1332 @*/ 1333 PetscErrorCode PCMGSetType(PC pc, PCMGType form) 1334 { 1335 PC_MG *mg = (PC_MG *)pc->data; 1336 1337 PetscFunctionBegin; 1338 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1339 PetscValidLogicalCollectiveEnum(pc, form, 2); 1340 mg->am = form; 1341 if (form == PC_MG_MULTIPLICATIVE) pc->ops->applyrichardson = PCApplyRichardson_MG; 1342 else pc->ops->applyrichardson = NULL; 1343 PetscFunctionReturn(PETSC_SUCCESS); 1344 } 1345 1346 /*@ 1347 PCMGGetType - Finds the form of multigrid the `PCMG` is using multiplicative, additive, full, or the Kaskade algorithm. 1348 1349 Logically Collective 1350 1351 Input Parameter: 1352 . pc - the preconditioner context 1353 1354 Output Parameter: 1355 . type - one of `PC_MG_MULTIPLICATIVE`, `PC_MG_ADDITIVE`, `PC_MG_FULL`, `PC_MG_KASKADE`, `PCMGCycleType` 1356 1357 Level: advanced 1358 1359 .seealso: [](ch_ksp), `PCMGType`, `PCMG`, `PCMGGetLevels()`, `PCMGSetLevels()`, `PCMGSetType()` 1360 @*/ 1361 PetscErrorCode PCMGGetType(PC pc, PCMGType *type) 1362 { 1363 PC_MG *mg = (PC_MG *)pc->data; 1364 1365 PetscFunctionBegin; 1366 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1367 *type = mg->am; 1368 PetscFunctionReturn(PETSC_SUCCESS); 1369 } 1370 1371 /*@ 1372 PCMGSetCycleType - Sets the type cycles to use. Use `PCMGSetCycleTypeOnLevel()` for more 1373 complicated cycling. 1374 1375 Logically Collective 1376 1377 Input Parameters: 1378 + pc - the multigrid context 1379 - n - either `PC_MG_CYCLE_V` or `PC_MG_CYCLE_W` 1380 1381 Options Database Key: 1382 . -pc_mg_cycle_type <v,w> - provide the cycle desired 1383 1384 Level: advanced 1385 1386 .seealso: [](ch_ksp), `PCMG`, `PCMGSetCycleTypeOnLevel()`, `PCMGType`, `PCMGCycleType` 1387 @*/ 1388 PetscErrorCode PCMGSetCycleType(PC pc, PCMGCycleType n) 1389 { 1390 PC_MG *mg = (PC_MG *)pc->data; 1391 PC_MG_Levels **mglevels = mg->levels; 1392 PetscInt i, levels; 1393 1394 PetscFunctionBegin; 1395 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1396 PetscValidLogicalCollectiveEnum(pc, n, 2); 1397 PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ORDER, "Must set MG levels with PCMGSetLevels() before calling"); 1398 levels = mglevels[0]->levels; 1399 for (i = 0; i < levels; i++) mglevels[i]->cycles = n; 1400 PetscFunctionReturn(PETSC_SUCCESS); 1401 } 1402 1403 /*@ 1404 PCMGMultiplicativeSetCycles - Sets the number of cycles to use for each preconditioner step 1405 of multigrid when `PCMGType` is `PC_MG_MULTIPLICATIVE` 1406 1407 Logically Collective 1408 1409 Input Parameters: 1410 + pc - the multigrid context 1411 - n - number of cycles (default is 1) 1412 1413 Options Database Key: 1414 . -pc_mg_multiplicative_cycles n - set the number of cycles 1415 1416 Level: advanced 1417 1418 Note: 1419 This is not associated with setting a v or w cycle, that is set with `PCMGSetCycleType()` 1420 1421 .seealso: [](ch_ksp), `PCMGSetCycleTypeOnLevel()`, `PCMGSetCycleType()`, `PCMGCycleType`, `PCMGType` 1422 @*/ 1423 PetscErrorCode PCMGMultiplicativeSetCycles(PC pc, PetscInt n) 1424 { 1425 PC_MG *mg = (PC_MG *)pc->data; 1426 1427 PetscFunctionBegin; 1428 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1429 PetscValidLogicalCollectiveInt(pc, n, 2); 1430 mg->cyclesperpcapply = n; 1431 PetscFunctionReturn(PETSC_SUCCESS); 1432 } 1433 1434 static PetscErrorCode PCMGSetGalerkin_MG(PC pc, PCMGGalerkinType use) 1435 { 1436 PC_MG *mg = (PC_MG *)pc->data; 1437 1438 PetscFunctionBegin; 1439 mg->galerkin = use; 1440 PetscFunctionReturn(PETSC_SUCCESS); 1441 } 1442 1443 /*@ 1444 PCMGSetGalerkin - Causes the coarser grid matrices to be computed from the 1445 finest grid via the Galerkin process: A_i-1 = r_i * A_i * p_i 1446 1447 Logically Collective 1448 1449 Input Parameters: 1450 + pc - the multigrid context 1451 - use - one of `PC_MG_GALERKIN_BOTH`, `PC_MG_GALERKIN_PMAT`, `PC_MG_GALERKIN_MAT`, or `PC_MG_GALERKIN_NONE` 1452 1453 Options Database Key: 1454 . -pc_mg_galerkin <both,pmat,mat,none> - set the matrices to form via the Galerkin process 1455 1456 Level: intermediate 1457 1458 Note: 1459 Some codes that use `PCMG` such as `PCGAMG` use Galerkin internally while constructing the hierarchy and thus do not 1460 use the `PCMG` construction of the coarser grids. 1461 1462 .seealso: [](ch_ksp), `PCMG`, `PCMGGetGalerkin()`, `PCMGGalerkinType` 1463 @*/ 1464 PetscErrorCode PCMGSetGalerkin(PC pc, PCMGGalerkinType use) 1465 { 1466 PetscFunctionBegin; 1467 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1468 PetscTryMethod(pc, "PCMGSetGalerkin_C", (PC, PCMGGalerkinType), (pc, use)); 1469 PetscFunctionReturn(PETSC_SUCCESS); 1470 } 1471 1472 /*@ 1473 PCMGGetGalerkin - Checks if Galerkin multigrid is being used, i.e. A_i-1 = r_i * A_i * p_i 1474 1475 Not Collective 1476 1477 Input Parameter: 1478 . pc - the multigrid context 1479 1480 Output Parameter: 1481 . galerkin - one of `PC_MG_GALERKIN_BOTH`,`PC_MG_GALERKIN_PMAT`,`PC_MG_GALERKIN_MAT`, `PC_MG_GALERKIN_NONE`, or `PC_MG_GALERKIN_EXTERNAL` 1482 1483 Level: intermediate 1484 1485 .seealso: [](ch_ksp), `PCMG`, `PCMGSetGalerkin()`, `PCMGGalerkinType` 1486 @*/ 1487 PetscErrorCode PCMGGetGalerkin(PC pc, PCMGGalerkinType *galerkin) 1488 { 1489 PC_MG *mg = (PC_MG *)pc->data; 1490 1491 PetscFunctionBegin; 1492 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1493 *galerkin = mg->galerkin; 1494 PetscFunctionReturn(PETSC_SUCCESS); 1495 } 1496 1497 static PetscErrorCode PCMGSetAdaptInterpolation_MG(PC pc, PetscBool adapt) 1498 { 1499 PC_MG *mg = (PC_MG *)pc->data; 1500 1501 PetscFunctionBegin; 1502 mg->adaptInterpolation = adapt; 1503 PetscFunctionReturn(PETSC_SUCCESS); 1504 } 1505 1506 static PetscErrorCode PCMGGetAdaptInterpolation_MG(PC pc, PetscBool *adapt) 1507 { 1508 PC_MG *mg = (PC_MG *)pc->data; 1509 1510 PetscFunctionBegin; 1511 *adapt = mg->adaptInterpolation; 1512 PetscFunctionReturn(PETSC_SUCCESS); 1513 } 1514 1515 static PetscErrorCode PCMGSetAdaptCoarseSpaceType_MG(PC pc, PCMGCoarseSpaceType ctype) 1516 { 1517 PC_MG *mg = (PC_MG *)pc->data; 1518 1519 PetscFunctionBegin; 1520 mg->adaptInterpolation = ctype != PCMG_ADAPT_NONE ? PETSC_TRUE : PETSC_FALSE; 1521 mg->coarseSpaceType = ctype; 1522 PetscCall(PCMGSetGalerkin(pc, PC_MG_GALERKIN_BOTH)); 1523 PetscFunctionReturn(PETSC_SUCCESS); 1524 } 1525 1526 static PetscErrorCode PCMGGetAdaptCoarseSpaceType_MG(PC pc, PCMGCoarseSpaceType *ctype) 1527 { 1528 PC_MG *mg = (PC_MG *)pc->data; 1529 1530 PetscFunctionBegin; 1531 *ctype = mg->coarseSpaceType; 1532 PetscFunctionReturn(PETSC_SUCCESS); 1533 } 1534 1535 static PetscErrorCode PCMGSetAdaptCR_MG(PC pc, PetscBool cr) 1536 { 1537 PC_MG *mg = (PC_MG *)pc->data; 1538 1539 PetscFunctionBegin; 1540 mg->compatibleRelaxation = cr; 1541 PetscFunctionReturn(PETSC_SUCCESS); 1542 } 1543 1544 static PetscErrorCode PCMGGetAdaptCR_MG(PC pc, PetscBool *cr) 1545 { 1546 PC_MG *mg = (PC_MG *)pc->data; 1547 1548 PetscFunctionBegin; 1549 *cr = mg->compatibleRelaxation; 1550 PetscFunctionReturn(PETSC_SUCCESS); 1551 } 1552 1553 /*@ 1554 PCMGSetAdaptCoarseSpaceType - Set the type of adaptive coarse space. 1555 1556 Adapts or creates the interpolator based upon a vector space which should be accurately captured by the next coarser mesh, and thus accurately interpolated. 1557 1558 Logically Collective 1559 1560 Input Parameters: 1561 + pc - the multigrid context 1562 - ctype - the type of coarse space 1563 1564 Options Database Keys: 1565 + -pc_mg_adapt_interp_n <int> - The number of modes to use 1566 - -pc_mg_adapt_interp_coarse_space <type> - The type of coarse space: none, `polynomial`, `harmonic`, `eigenvector`, `generalized_eigenvector`, `gdsw` 1567 1568 Level: intermediate 1569 1570 Note: 1571 Requires a `DM` with specific functionality be attached to the `PC`. 1572 1573 .seealso: [](ch_ksp), `PCMG`, `PCMGCoarseSpaceType`, `PCMGGetAdaptCoarseSpaceType()`, `PCMGSetGalerkin()`, `PCMGSetAdaptInterpolation()`, `DM` 1574 @*/ 1575 PetscErrorCode PCMGSetAdaptCoarseSpaceType(PC pc, PCMGCoarseSpaceType ctype) 1576 { 1577 PetscFunctionBegin; 1578 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1579 PetscValidLogicalCollectiveEnum(pc, ctype, 2); 1580 PetscTryMethod(pc, "PCMGSetAdaptCoarseSpaceType_C", (PC, PCMGCoarseSpaceType), (pc, ctype)); 1581 PetscFunctionReturn(PETSC_SUCCESS); 1582 } 1583 1584 /*@ 1585 PCMGGetAdaptCoarseSpaceType - Get the type of adaptive coarse space. 1586 1587 Not Collective 1588 1589 Input Parameter: 1590 . pc - the multigrid context 1591 1592 Output Parameter: 1593 . ctype - the type of coarse space 1594 1595 Level: intermediate 1596 1597 .seealso: [](ch_ksp), `PCMG`, `PCMGCoarseSpaceType`, `PCMGSetAdaptCoarseSpaceType()`, `PCMGSetGalerkin()`, `PCMGSetAdaptInterpolation()` 1598 @*/ 1599 PetscErrorCode PCMGGetAdaptCoarseSpaceType(PC pc, PCMGCoarseSpaceType *ctype) 1600 { 1601 PetscFunctionBegin; 1602 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1603 PetscAssertPointer(ctype, 2); 1604 PetscUseMethod(pc, "PCMGGetAdaptCoarseSpaceType_C", (PC, PCMGCoarseSpaceType *), (pc, ctype)); 1605 PetscFunctionReturn(PETSC_SUCCESS); 1606 } 1607 1608 /* MATT: REMOVE? */ 1609 /*@ 1610 PCMGSetAdaptInterpolation - Adapt the interpolator based upon a vector space which should be accurately captured by the next coarser mesh, and thus accurately interpolated. 1611 1612 Logically Collective 1613 1614 Input Parameters: 1615 + pc - the multigrid context 1616 - adapt - flag for adaptation of the interpolator 1617 1618 Options Database Keys: 1619 + -pc_mg_adapt_interp - Turn on adaptation 1620 . -pc_mg_adapt_interp_n <int> - The number of modes to use, should be divisible by dimension 1621 - -pc_mg_adapt_interp_coarse_space <type> - The type of coarse space: polynomial, harmonic, eigenvector, generalized_eigenvector 1622 1623 Level: intermediate 1624 1625 .seealso: [](ch_ksp), `PCMG`, `PCMGGetAdaptInterpolation()`, `PCMGSetGalerkin()`, `PCMGGetAdaptCoarseSpaceType()`, `PCMGSetAdaptCoarseSpaceType()` 1626 @*/ 1627 PetscErrorCode PCMGSetAdaptInterpolation(PC pc, PetscBool adapt) 1628 { 1629 PetscFunctionBegin; 1630 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1631 PetscTryMethod(pc, "PCMGSetAdaptInterpolation_C", (PC, PetscBool), (pc, adapt)); 1632 PetscFunctionReturn(PETSC_SUCCESS); 1633 } 1634 1635 /*@ 1636 PCMGGetAdaptInterpolation - Get the flag to adapt the interpolator based upon a vector space which should be accurately captured by the next coarser mesh, 1637 and thus accurately interpolated. 1638 1639 Not Collective 1640 1641 Input Parameter: 1642 . pc - the multigrid context 1643 1644 Output Parameter: 1645 . adapt - flag for adaptation of the interpolator 1646 1647 Level: intermediate 1648 1649 .seealso: [](ch_ksp), `PCMG`, `PCMGSetAdaptInterpolation()`, `PCMGSetGalerkin()`, `PCMGGetAdaptCoarseSpaceType()`, `PCMGSetAdaptCoarseSpaceType()` 1650 @*/ 1651 PetscErrorCode PCMGGetAdaptInterpolation(PC pc, PetscBool *adapt) 1652 { 1653 PetscFunctionBegin; 1654 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1655 PetscAssertPointer(adapt, 2); 1656 PetscUseMethod(pc, "PCMGGetAdaptInterpolation_C", (PC, PetscBool *), (pc, adapt)); 1657 PetscFunctionReturn(PETSC_SUCCESS); 1658 } 1659 1660 /*@ 1661 PCMGSetAdaptCR - Monitor the coarse space quality using an auxiliary solve with compatible relaxation. 1662 1663 Logically Collective 1664 1665 Input Parameters: 1666 + pc - the multigrid context 1667 - cr - flag for compatible relaxation 1668 1669 Options Database Key: 1670 . -pc_mg_adapt_cr - Turn on compatible relaxation 1671 1672 Level: intermediate 1673 1674 .seealso: [](ch_ksp), `PCMG`, `PCMGGetAdaptCR()`, `PCMGSetAdaptInterpolation()`, `PCMGSetGalerkin()`, `PCMGGetAdaptCoarseSpaceType()`, `PCMGSetAdaptCoarseSpaceType()` 1675 @*/ 1676 PetscErrorCode PCMGSetAdaptCR(PC pc, PetscBool cr) 1677 { 1678 PetscFunctionBegin; 1679 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1680 PetscTryMethod(pc, "PCMGSetAdaptCR_C", (PC, PetscBool), (pc, cr)); 1681 PetscFunctionReturn(PETSC_SUCCESS); 1682 } 1683 1684 /*@ 1685 PCMGGetAdaptCR - Get the flag to monitor coarse space quality using an auxiliary solve with compatible relaxation. 1686 1687 Not Collective 1688 1689 Input Parameter: 1690 . pc - the multigrid context 1691 1692 Output Parameter: 1693 . cr - flag for compatible relaxaion 1694 1695 Level: intermediate 1696 1697 .seealso: [](ch_ksp), `PCMGSetAdaptCR()`, `PCMGGetAdaptInterpolation()`, `PCMGSetGalerkin()`, `PCMGGetAdaptCoarseSpaceType()`, `PCMGSetAdaptCoarseSpaceType()` 1698 @*/ 1699 PetscErrorCode PCMGGetAdaptCR(PC pc, PetscBool *cr) 1700 { 1701 PetscFunctionBegin; 1702 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1703 PetscAssertPointer(cr, 2); 1704 PetscUseMethod(pc, "PCMGGetAdaptCR_C", (PC, PetscBool *), (pc, cr)); 1705 PetscFunctionReturn(PETSC_SUCCESS); 1706 } 1707 1708 /*@ 1709 PCMGSetNumberSmooth - Sets the number of pre and post-smoothing steps to use 1710 on all levels. Use `PCMGDistinctSmoothUp()` to create separate up and down smoothers if you want different numbers of 1711 pre- and post-smoothing steps. 1712 1713 Logically Collective 1714 1715 Input Parameters: 1716 + pc - the multigrid context 1717 - n - the number of smoothing steps 1718 1719 Options Database Key: 1720 . -mg_levels_ksp_max_it <n> - Sets number of pre and post-smoothing steps 1721 1722 Level: advanced 1723 1724 Note: 1725 This does not set a value on the coarsest grid, since we assume that there is no separate smooth up on the coarsest grid. 1726 1727 .seealso: [](ch_ksp), `PCMG`, `PCMGSetDistinctSmoothUp()` 1728 @*/ 1729 PetscErrorCode PCMGSetNumberSmooth(PC pc, PetscInt n) 1730 { 1731 PC_MG *mg = (PC_MG *)pc->data; 1732 PC_MG_Levels **mglevels = mg->levels; 1733 PetscInt i, levels; 1734 1735 PetscFunctionBegin; 1736 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1737 PetscValidLogicalCollectiveInt(pc, n, 2); 1738 PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ORDER, "Must set MG levels with PCMGSetLevels() before calling"); 1739 levels = mglevels[0]->levels; 1740 1741 for (i = 1; i < levels; i++) { 1742 PetscCall(KSPSetTolerances(mglevels[i]->smoothu, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT, n)); 1743 PetscCall(KSPSetTolerances(mglevels[i]->smoothd, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT, n)); 1744 mg->default_smoothu = n; 1745 mg->default_smoothd = n; 1746 } 1747 PetscFunctionReturn(PETSC_SUCCESS); 1748 } 1749 1750 /*@ 1751 PCMGSetDistinctSmoothUp - sets the up (post) smoother to be a separate `KSP` from the down (pre) smoother on all levels 1752 and adds the suffix _up to the options name 1753 1754 Logically Collective 1755 1756 Input Parameter: 1757 . pc - the preconditioner context 1758 1759 Options Database Key: 1760 . -pc_mg_distinct_smoothup <bool> - use distinct smoothing objects 1761 1762 Level: advanced 1763 1764 Note: 1765 This does not set a value on the coarsest grid, since we assume that there is no separate smooth up on the coarsest grid. 1766 1767 .seealso: [](ch_ksp), `PCMG`, `PCMGSetNumberSmooth()` 1768 @*/ 1769 PetscErrorCode PCMGSetDistinctSmoothUp(PC pc) 1770 { 1771 PC_MG *mg = (PC_MG *)pc->data; 1772 PC_MG_Levels **mglevels = mg->levels; 1773 PetscInt i, levels; 1774 KSP subksp; 1775 1776 PetscFunctionBegin; 1777 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1778 PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ORDER, "Must set MG levels with PCMGSetLevels() before calling"); 1779 levels = mglevels[0]->levels; 1780 1781 for (i = 1; i < levels; i++) { 1782 const char *prefix = NULL; 1783 /* make sure smoother up and down are different */ 1784 PetscCall(PCMGGetSmootherUp(pc, i, &subksp)); 1785 PetscCall(KSPGetOptionsPrefix(mglevels[i]->smoothd, &prefix)); 1786 PetscCall(KSPSetOptionsPrefix(subksp, prefix)); 1787 PetscCall(KSPAppendOptionsPrefix(subksp, "up_")); 1788 } 1789 PetscFunctionReturn(PETSC_SUCCESS); 1790 } 1791 1792 /* No new matrices are created, and the coarse operator matrices are the references to the original ones */ 1793 static PetscErrorCode PCGetInterpolations_MG(PC pc, PetscInt *num_levels, Mat *interpolations[]) 1794 { 1795 PC_MG *mg = (PC_MG *)pc->data; 1796 PC_MG_Levels **mglevels = mg->levels; 1797 Mat *mat; 1798 PetscInt l; 1799 1800 PetscFunctionBegin; 1801 PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set MG levels before calling"); 1802 PetscCall(PetscMalloc1(mg->nlevels, &mat)); 1803 for (l = 1; l < mg->nlevels; l++) { 1804 mat[l - 1] = mglevels[l]->interpolate; 1805 PetscCall(PetscObjectReference((PetscObject)mat[l - 1])); 1806 } 1807 *num_levels = mg->nlevels; 1808 *interpolations = mat; 1809 PetscFunctionReturn(PETSC_SUCCESS); 1810 } 1811 1812 /* No new matrices are created, and the coarse operator matrices are the references to the original ones */ 1813 static PetscErrorCode PCGetCoarseOperators_MG(PC pc, PetscInt *num_levels, Mat *coarseOperators[]) 1814 { 1815 PC_MG *mg = (PC_MG *)pc->data; 1816 PC_MG_Levels **mglevels = mg->levels; 1817 PetscInt l; 1818 Mat *mat; 1819 1820 PetscFunctionBegin; 1821 PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set MG levels before calling"); 1822 PetscCall(PetscMalloc1(mg->nlevels, &mat)); 1823 for (l = 0; l < mg->nlevels - 1; l++) { 1824 PetscCall(KSPGetOperators(mglevels[l]->smoothd, NULL, &mat[l])); 1825 PetscCall(PetscObjectReference((PetscObject)mat[l])); 1826 } 1827 *num_levels = mg->nlevels; 1828 *coarseOperators = mat; 1829 PetscFunctionReturn(PETSC_SUCCESS); 1830 } 1831 1832 /*@C 1833 PCMGRegisterCoarseSpaceConstructor - Adds a method to the `PCMG` package for coarse space construction. 1834 1835 Not Collective, No Fortran Support 1836 1837 Input Parameters: 1838 + name - name of the constructor 1839 - function - constructor routine 1840 1841 Calling sequence of `function`: 1842 + pc - The `PC` object 1843 . l - The multigrid level, 0 is the coarse level 1844 . dm - The `DM` for this level 1845 . smooth - The level smoother 1846 . Nc - The size of the coarse space 1847 . initGuess - Basis for an initial guess for the space 1848 - coarseSp - A basis for the computed coarse space 1849 1850 Level: advanced 1851 1852 Developer Notes: 1853 How come this is not used by `PCGAMG`? 1854 1855 .seealso: [](ch_ksp), `PCMG`, `PCMGGetCoarseSpaceConstructor()`, `PCRegister()` 1856 @*/ 1857 PetscErrorCode PCMGRegisterCoarseSpaceConstructor(const char name[], PetscErrorCode (*function)(PC pc, PetscInt l, DM dm, KSP smooth, PetscInt Nc, Mat initGuess, Mat *coarseSp)) 1858 { 1859 PetscFunctionBegin; 1860 PetscCall(PCInitializePackage()); 1861 PetscCall(PetscFunctionListAdd(&PCMGCoarseList, name, function)); 1862 PetscFunctionReturn(PETSC_SUCCESS); 1863 } 1864 1865 /*@C 1866 PCMGGetCoarseSpaceConstructor - Returns the given coarse space construction method. 1867 1868 Not Collective, No Fortran Support 1869 1870 Input Parameter: 1871 . name - name of the constructor 1872 1873 Output Parameter: 1874 . function - constructor routine 1875 1876 Level: advanced 1877 1878 .seealso: [](ch_ksp), `PCMG`, `PCMGRegisterCoarseSpaceConstructor()`, `PCRegister()` 1879 @*/ 1880 PetscErrorCode PCMGGetCoarseSpaceConstructor(const char name[], PetscErrorCode (**function)(PC, PetscInt, DM, KSP, PetscInt, Mat, Mat *)) 1881 { 1882 PetscFunctionBegin; 1883 PetscCall(PetscFunctionListFind(PCMGCoarseList, name, function)); 1884 PetscFunctionReturn(PETSC_SUCCESS); 1885 } 1886 1887 /*MC 1888 PCMG - Use multigrid preconditioning. This preconditioner requires you provide additional 1889 information about the coarser grid matrices and restriction/interpolation operators. 1890 1891 Options Database Keys: 1892 + -pc_mg_levels <nlevels> - number of levels including finest 1893 . -pc_mg_cycle_type <v,w> - provide the cycle desired 1894 . -pc_mg_type <additive,multiplicative,full,kaskade> - multiplicative is the default 1895 . -pc_mg_log - log information about time spent on each level of the solver 1896 . -pc_mg_distinct_smoothup - configure up (after interpolation) and down (before restriction) smoothers separately (with different options prefixes) 1897 . -pc_mg_galerkin <both,pmat,mat,none> - use Galerkin process to compute coarser operators, i.e. Acoarse = R A R' 1898 . -pc_mg_multiplicative_cycles - number of cycles to use as the preconditioner (defaults to 1) 1899 . -pc_mg_dump_matlab - dumps the matrices for each level and the restriction/interpolation matrices 1900 to a `PETSCVIEWERSOCKET` for reading from MATLAB. 1901 - -pc_mg_dump_binary -dumps the matrices for each level and the restriction/interpolation matrices 1902 to the binary output file called binaryoutput 1903 1904 Level: intermediate 1905 1906 Notes: 1907 The Krylov solver (if any) and preconditioner (smoother) and their parameters are controlled from the options database with the standard 1908 options database keywords prefixed with `-mg_levels_` to affect all the levels but the coarsest, which is controlled with `-mg_coarse_`, 1909 and the finest where `-mg_fine_` can override `-mg_levels_`. One can set different preconditioners etc on specific levels with the prefix 1910 `-mg_levels_n_` where `n` is the level number (zero being the coarse level. For example 1911 .vb 1912 -mg_levels_ksp_type gmres -mg_levels_pc_type bjacobi -mg_coarse_pc_type svd -mg_levels_2_pc_type sor 1913 .ve 1914 These options also work for controlling the smoothers etc inside `PCGAMG` 1915 1916 If one uses a Krylov method such `KSPGMRES` or `KSPCG` as the smoother then one must use `KSPFGMRES`, `KSPGCR`, or `KSPRICHARDSON` as the outer Krylov method 1917 1918 When run with a single level the smoother options are used on that level NOT the coarse grid solver options 1919 1920 When run with `KSPRICHARDSON` the convergence test changes slightly if monitor is turned on. The iteration count may change slightly. This 1921 is because without monitoring the residual norm is computed WITHIN each multigrid cycle on the finest level after the pre-smoothing 1922 (because the residual has just been computed for the multigrid algorithm and is hence available for free) while with monitoring the 1923 residual is computed at the end of each cycle. 1924 1925 .seealso: [](sec_mg), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCMGType`, `PCEXOTIC`, `PCGAMG`, `PCML`, `PCHYPRE` 1926 `PCMGSetLevels()`, `PCMGGetLevels()`, `PCMGSetType()`, `PCMGSetCycleType()`, 1927 `PCMGSetDistinctSmoothUp()`, `PCMGGetCoarseSolve()`, `PCMGSetResidual()`, `PCMGSetInterpolation()`, 1928 `PCMGSetRestriction()`, `PCMGGetSmoother()`, `PCMGGetSmootherUp()`, `PCMGGetSmootherDown()`, 1929 `PCMGSetCycleTypeOnLevel()`, `PCMGSetRhs()`, `PCMGSetX()`, `PCMGSetR()`, 1930 `PCMGSetAdaptCR()`, `PCMGGetAdaptInterpolation()`, `PCMGSetGalerkin()`, `PCMGGetAdaptCoarseSpaceType()`, `PCMGSetAdaptCoarseSpaceType()` 1931 M*/ 1932 1933 PETSC_EXTERN PetscErrorCode PCCreate_MG(PC pc) 1934 { 1935 PC_MG *mg; 1936 1937 PetscFunctionBegin; 1938 PetscCall(PetscNew(&mg)); 1939 pc->data = mg; 1940 mg->nlevels = -1; 1941 mg->am = PC_MG_MULTIPLICATIVE; 1942 mg->galerkin = PC_MG_GALERKIN_NONE; 1943 mg->adaptInterpolation = PETSC_FALSE; 1944 mg->Nc = -1; 1945 mg->eigenvalue = -1; 1946 1947 pc->useAmat = PETSC_TRUE; 1948 1949 pc->ops->apply = PCApply_MG; 1950 pc->ops->applytranspose = PCApplyTranspose_MG; 1951 pc->ops->matapply = PCMatApply_MG; 1952 pc->ops->setup = PCSetUp_MG; 1953 pc->ops->reset = PCReset_MG; 1954 pc->ops->destroy = PCDestroy_MG; 1955 pc->ops->setfromoptions = PCSetFromOptions_MG; 1956 pc->ops->view = PCView_MG; 1957 1958 PetscCall(PetscObjectComposedDataRegister(&mg->eigenvalue)); 1959 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGSetGalerkin_C", PCMGSetGalerkin_MG)); 1960 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGetLevels_C", PCMGGetLevels_MG)); 1961 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGSetLevels_C", PCMGSetLevels_MG)); 1962 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGetInterpolations_C", PCGetInterpolations_MG)); 1963 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGetCoarseOperators_C", PCGetCoarseOperators_MG)); 1964 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGSetAdaptInterpolation_C", PCMGSetAdaptInterpolation_MG)); 1965 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGetAdaptInterpolation_C", PCMGGetAdaptInterpolation_MG)); 1966 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGSetAdaptCR_C", PCMGSetAdaptCR_MG)); 1967 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGetAdaptCR_C", PCMGGetAdaptCR_MG)); 1968 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGSetAdaptCoarseSpaceType_C", PCMGSetAdaptCoarseSpaceType_MG)); 1969 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGetAdaptCoarseSpaceType_C", PCMGGetAdaptCoarseSpaceType_MG)); 1970 PetscFunctionReturn(PETSC_SUCCESS); 1971 } 1972