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