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 Parameter: 1334 + gc - operator complexity = sum_i(nnz_i) / nnz_0 1335 . oc - grid complexity = sum_i(n_i) / n_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 if (!pc->setupcalled) { 1352 *gc = *oc = 0; 1353 PetscFunctionReturn(0); 1354 } 1355 PetscCheck(mg->nlevels > 0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"MG has no levels"); 1356 for (lev=0; lev<mg->nlevels; lev++) { 1357 Mat dB; 1358 ierr = KSPGetOperators(mglevels[lev]->smoothd,NULL,&dB);CHKERRQ(ierr); 1359 ierr = MatGetInfo(dB,MAT_GLOBAL_SUM,&info);CHKERRQ(ierr); /* global reduction */ 1360 ierr = MatGetSize(dB,&N,NULL);CHKERRQ(ierr); 1361 sgc += N; 1362 soc += info.nz_used; 1363 if (lev==mg->nlevels-1) {nnz0 = info.nz_used; n0 = N;} 1364 } 1365 if (n0 > 0) *gc = (PetscReal)(sgc/n0); 1366 else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Number for grid points on finest level is not available"); 1367 if (nnz0 > 0) *oc = (PetscReal)(soc/nnz0); 1368 PetscFunctionReturn(0); 1369 } 1370 1371 /*@ 1372 PCMGSetType - Determines the form of multigrid to use: 1373 multiplicative, additive, full, or the Kaskade algorithm. 1374 1375 Logically Collective on PC 1376 1377 Input Parameters: 1378 + pc - the preconditioner context 1379 - form - multigrid form, one of PC_MG_MULTIPLICATIVE, PC_MG_ADDITIVE, 1380 PC_MG_FULL, PC_MG_KASKADE 1381 1382 Options Database Key: 1383 . -pc_mg_type <form> - Sets <form>, one of multiplicative, 1384 additive, full, kaskade 1385 1386 Level: advanced 1387 1388 .seealso: PCMGSetLevels() 1389 @*/ 1390 PetscErrorCode PCMGSetType(PC pc,PCMGType form) 1391 { 1392 PC_MG *mg = (PC_MG*)pc->data; 1393 1394 PetscFunctionBegin; 1395 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1396 PetscValidLogicalCollectiveEnum(pc,form,2); 1397 mg->am = form; 1398 if (form == PC_MG_MULTIPLICATIVE) pc->ops->applyrichardson = PCApplyRichardson_MG; 1399 else pc->ops->applyrichardson = NULL; 1400 PetscFunctionReturn(0); 1401 } 1402 1403 /*@ 1404 PCMGGetType - Determines the form of multigrid to use: 1405 multiplicative, additive, full, or the Kaskade algorithm. 1406 1407 Logically Collective on PC 1408 1409 Input Parameter: 1410 . pc - the preconditioner context 1411 1412 Output Parameter: 1413 . type - one of PC_MG_MULTIPLICATIVE, PC_MG_ADDITIVE,PC_MG_FULL, PC_MG_KASKADE 1414 1415 Level: advanced 1416 1417 .seealso: PCMGSetLevels() 1418 @*/ 1419 PetscErrorCode PCMGGetType(PC pc,PCMGType *type) 1420 { 1421 PC_MG *mg = (PC_MG*)pc->data; 1422 1423 PetscFunctionBegin; 1424 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1425 *type = mg->am; 1426 PetscFunctionReturn(0); 1427 } 1428 1429 /*@ 1430 PCMGSetCycleType - Sets the type cycles to use. Use PCMGSetCycleTypeOnLevel() for more 1431 complicated cycling. 1432 1433 Logically Collective on PC 1434 1435 Input Parameters: 1436 + pc - the multigrid context 1437 - n - either PC_MG_CYCLE_V or PC_MG_CYCLE_W 1438 1439 Options Database Key: 1440 . -pc_mg_cycle_type <v,w> - provide the cycle desired 1441 1442 Level: advanced 1443 1444 .seealso: PCMGSetCycleTypeOnLevel() 1445 @*/ 1446 PetscErrorCode PCMGSetCycleType(PC pc,PCMGCycleType n) 1447 { 1448 PC_MG *mg = (PC_MG*)pc->data; 1449 PC_MG_Levels **mglevels = mg->levels; 1450 PetscInt i,levels; 1451 1452 PetscFunctionBegin; 1453 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1454 PetscValidLogicalCollectiveEnum(pc,n,2); 1455 PetscCheckFalse(!mglevels,PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Must set MG levels with PCMGSetLevels() before calling"); 1456 levels = mglevels[0]->levels; 1457 for (i=0; i<levels; i++) mglevels[i]->cycles = n; 1458 PetscFunctionReturn(0); 1459 } 1460 1461 /*@ 1462 PCMGMultiplicativeSetCycles - Sets the number of cycles to use for each preconditioner step 1463 of multigrid when PCMGType of PC_MG_MULTIPLICATIVE is used 1464 1465 Logically Collective on PC 1466 1467 Input Parameters: 1468 + pc - the multigrid context 1469 - n - number of cycles (default is 1) 1470 1471 Options Database Key: 1472 . -pc_mg_multiplicative_cycles n 1473 1474 Level: advanced 1475 1476 Notes: 1477 This is not associated with setting a v or w cycle, that is set with PCMGSetCycleType() 1478 1479 .seealso: PCMGSetCycleTypeOnLevel(), PCMGSetCycleType() 1480 @*/ 1481 PetscErrorCode PCMGMultiplicativeSetCycles(PC pc,PetscInt n) 1482 { 1483 PC_MG *mg = (PC_MG*)pc->data; 1484 1485 PetscFunctionBegin; 1486 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1487 PetscValidLogicalCollectiveInt(pc,n,2); 1488 mg->cyclesperpcapply = n; 1489 PetscFunctionReturn(0); 1490 } 1491 1492 PetscErrorCode PCMGSetGalerkin_MG(PC pc,PCMGGalerkinType use) 1493 { 1494 PC_MG *mg = (PC_MG*)pc->data; 1495 1496 PetscFunctionBegin; 1497 mg->galerkin = use; 1498 PetscFunctionReturn(0); 1499 } 1500 1501 /*@ 1502 PCMGSetGalerkin - Causes the coarser grid matrices to be computed from the 1503 finest grid via the Galerkin process: A_i-1 = r_i * A_i * p_i 1504 1505 Logically Collective on PC 1506 1507 Input Parameters: 1508 + pc - the multigrid context 1509 - use - one of PC_MG_GALERKIN_BOTH,PC_MG_GALERKIN_PMAT,PC_MG_GALERKIN_MAT, or PC_MG_GALERKIN_NONE 1510 1511 Options Database Key: 1512 . -pc_mg_galerkin <both,pmat,mat,none> 1513 1514 Level: intermediate 1515 1516 Notes: 1517 Some codes that use PCMG such as PCGAMG use Galerkin internally while constructing the hierarchy and thus do not 1518 use the PCMG construction of the coarser grids. 1519 1520 .seealso: PCMGGetGalerkin(), PCMGGalerkinType 1521 1522 @*/ 1523 PetscErrorCode PCMGSetGalerkin(PC pc,PCMGGalerkinType use) 1524 { 1525 PetscErrorCode ierr; 1526 1527 PetscFunctionBegin; 1528 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1529 ierr = PetscTryMethod(pc,"PCMGSetGalerkin_C",(PC,PCMGGalerkinType),(pc,use));CHKERRQ(ierr); 1530 PetscFunctionReturn(0); 1531 } 1532 1533 /*@ 1534 PCMGGetGalerkin - Checks if Galerkin multigrid is being used, i.e. 1535 A_i-1 = r_i * A_i * p_i 1536 1537 Not Collective 1538 1539 Input Parameter: 1540 . pc - the multigrid context 1541 1542 Output Parameter: 1543 . galerkin - one of PC_MG_GALERKIN_BOTH,PC_MG_GALERKIN_PMAT,PC_MG_GALERKIN_MAT, PC_MG_GALERKIN_NONE, or PC_MG_GALERKIN_EXTERNAL 1544 1545 Level: intermediate 1546 1547 .seealso: PCMGSetGalerkin(), PCMGGalerkinType 1548 1549 @*/ 1550 PetscErrorCode PCMGGetGalerkin(PC pc,PCMGGalerkinType *galerkin) 1551 { 1552 PC_MG *mg = (PC_MG*)pc->data; 1553 1554 PetscFunctionBegin; 1555 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1556 *galerkin = mg->galerkin; 1557 PetscFunctionReturn(0); 1558 } 1559 1560 PetscErrorCode PCMGSetAdaptInterpolation_MG(PC pc, PetscBool adapt) 1561 { 1562 PC_MG *mg = (PC_MG *) pc->data; 1563 1564 PetscFunctionBegin; 1565 mg->adaptInterpolation = adapt; 1566 PetscFunctionReturn(0); 1567 } 1568 1569 PetscErrorCode PCMGGetAdaptInterpolation_MG(PC pc, PetscBool *adapt) 1570 { 1571 PC_MG *mg = (PC_MG *) pc->data; 1572 1573 PetscFunctionBegin; 1574 *adapt = mg->adaptInterpolation; 1575 PetscFunctionReturn(0); 1576 } 1577 1578 PetscErrorCode PCMGSetAdaptCR_MG(PC pc, PetscBool cr) 1579 { 1580 PC_MG *mg = (PC_MG *) pc->data; 1581 1582 PetscFunctionBegin; 1583 mg->compatibleRelaxation = cr; 1584 PetscFunctionReturn(0); 1585 } 1586 1587 PetscErrorCode PCMGGetAdaptCR_MG(PC pc, PetscBool *cr) 1588 { 1589 PC_MG *mg = (PC_MG *) pc->data; 1590 1591 PetscFunctionBegin; 1592 *cr = mg->compatibleRelaxation; 1593 PetscFunctionReturn(0); 1594 } 1595 1596 /*@ 1597 PCMGSetAdaptInterpolation - Adapt the interpolator based upon a vector space which should be accurately captured by the next coarser mesh, and thus accurately interpolated. 1598 1599 Logically Collective on PC 1600 1601 Input Parameters: 1602 + pc - the multigrid context 1603 - adapt - flag for adaptation of the interpolator 1604 1605 Options Database Keys: 1606 + -pc_mg_adapt_interp - Turn on adaptation 1607 . -pc_mg_adapt_interp_n <int> - The number of modes to use, should be divisible by dimension 1608 - -pc_mg_adapt_interp_coarse_space <type> - The type of coarse space: polynomial, harmonic, eigenvector, generalized_eigenvector 1609 1610 Level: intermediate 1611 1612 .keywords: MG, set, Galerkin 1613 .seealso: PCMGGetAdaptInterpolation(), PCMGSetGalerkin() 1614 @*/ 1615 PetscErrorCode PCMGSetAdaptInterpolation(PC pc, PetscBool adapt) 1616 { 1617 PetscErrorCode ierr; 1618 1619 PetscFunctionBegin; 1620 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1621 ierr = PetscTryMethod(pc,"PCMGSetAdaptInterpolation_C",(PC,PetscBool),(pc,adapt));CHKERRQ(ierr); 1622 PetscFunctionReturn(0); 1623 } 1624 1625 /*@ 1626 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. 1627 1628 Not collective 1629 1630 Input Parameter: 1631 . pc - the multigrid context 1632 1633 Output Parameter: 1634 . adapt - flag for adaptation of the interpolator 1635 1636 Level: intermediate 1637 1638 .keywords: MG, set, Galerkin 1639 .seealso: PCMGSetAdaptInterpolation(), PCMGSetGalerkin() 1640 @*/ 1641 PetscErrorCode PCMGGetAdaptInterpolation(PC pc, PetscBool *adapt) 1642 { 1643 PetscErrorCode ierr; 1644 1645 PetscFunctionBegin; 1646 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1647 PetscValidPointer(adapt, 2); 1648 ierr = PetscUseMethod(pc,"PCMGGetAdaptInterpolation_C",(PC,PetscBool*),(pc,adapt));CHKERRQ(ierr); 1649 PetscFunctionReturn(0); 1650 } 1651 1652 /*@ 1653 PCMGSetAdaptCR - Monitor the coarse space quality using an auxiliary solve with compatible relaxation. 1654 1655 Logically Collective on PC 1656 1657 Input Parameters: 1658 + pc - the multigrid context 1659 - cr - flag for compatible relaxation 1660 1661 Options Database Keys: 1662 . -pc_mg_adapt_cr - Turn on compatible relaxation 1663 1664 Level: intermediate 1665 1666 .keywords: MG, set, Galerkin 1667 .seealso: PCMGGetAdaptCR(), PCMGSetAdaptInterpolation(), PCMGSetGalerkin() 1668 @*/ 1669 PetscErrorCode PCMGSetAdaptCR(PC pc, PetscBool cr) 1670 { 1671 PetscErrorCode ierr; 1672 1673 PetscFunctionBegin; 1674 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1675 ierr = PetscTryMethod(pc,"PCMGSetAdaptCR_C",(PC,PetscBool),(pc,cr));CHKERRQ(ierr); 1676 PetscFunctionReturn(0); 1677 } 1678 1679 /*@ 1680 PCMGGetAdaptCR - Get the flag to monitor coarse space quality using an auxiliary solve with compatible relaxation. 1681 1682 Not collective 1683 1684 Input Parameter: 1685 . pc - the multigrid context 1686 1687 Output Parameter: 1688 . cr - flag for compatible relaxaion 1689 1690 Level: intermediate 1691 1692 .keywords: MG, set, Galerkin 1693 .seealso: PCMGSetAdaptCR(), PCMGGetAdaptInterpolation(), PCMGSetGalerkin() 1694 @*/ 1695 PetscErrorCode PCMGGetAdaptCR(PC pc, PetscBool *cr) 1696 { 1697 PetscErrorCode ierr; 1698 1699 PetscFunctionBegin; 1700 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1701 PetscValidPointer(cr, 2); 1702 ierr = PetscUseMethod(pc,"PCMGGetAdaptCR_C",(PC,PetscBool*),(pc,cr));CHKERRQ(ierr); 1703 PetscFunctionReturn(0); 1704 } 1705 1706 /*@ 1707 PCMGSetNumberSmooth - Sets the number of pre and post-smoothing steps to use 1708 on all levels. Use PCMGDistinctSmoothUp() to create separate up and down smoothers if you want different numbers of 1709 pre- and post-smoothing steps. 1710 1711 Logically Collective on PC 1712 1713 Input Parameters: 1714 + mg - the multigrid context 1715 - n - the number of smoothing steps 1716 1717 Options Database Key: 1718 . -mg_levels_ksp_max_it <n> - Sets number of pre and post-smoothing steps 1719 1720 Level: advanced 1721 1722 Notes: 1723 this does not set a value on the coarsest grid, since we assume that 1724 there is no separate smooth up on the coarsest grid. 1725 1726 .seealso: PCMGSetDistinctSmoothUp() 1727 @*/ 1728 PetscErrorCode PCMGSetNumberSmooth(PC pc,PetscInt n) 1729 { 1730 PC_MG *mg = (PC_MG*)pc->data; 1731 PC_MG_Levels **mglevels = mg->levels; 1732 PetscErrorCode ierr; 1733 PetscInt i,levels; 1734 1735 PetscFunctionBegin; 1736 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1737 PetscValidLogicalCollectiveInt(pc,n,2); 1738 PetscCheckFalse(!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 ierr = KSPSetTolerances(mglevels[i]->smoothu,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,n);CHKERRQ(ierr); 1743 ierr = KSPSetTolerances(mglevels[i]->smoothd,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,n);CHKERRQ(ierr); 1744 mg->default_smoothu = n; 1745 mg->default_smoothd = n; 1746 } 1747 PetscFunctionReturn(0); 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 on PC 1755 1756 Input Parameters: 1757 . pc - the preconditioner context 1758 1759 Options Database Key: 1760 . -pc_mg_distinct_smoothup 1761 1762 Level: advanced 1763 1764 Notes: 1765 this does not set a value on the coarsest grid, since we assume that 1766 there is no separate smooth up on the coarsest grid. 1767 1768 .seealso: PCMGSetNumberSmooth() 1769 @*/ 1770 PetscErrorCode PCMGSetDistinctSmoothUp(PC pc) 1771 { 1772 PC_MG *mg = (PC_MG*)pc->data; 1773 PC_MG_Levels **mglevels = mg->levels; 1774 PetscErrorCode ierr; 1775 PetscInt i,levels; 1776 KSP subksp; 1777 1778 PetscFunctionBegin; 1779 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1780 PetscCheckFalse(!mglevels,PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Must set MG levels with PCMGSetLevels() before calling"); 1781 levels = mglevels[0]->levels; 1782 1783 for (i=1; i<levels; i++) { 1784 const char *prefix = NULL; 1785 /* make sure smoother up and down are different */ 1786 ierr = PCMGGetSmootherUp(pc,i,&subksp);CHKERRQ(ierr); 1787 ierr = KSPGetOptionsPrefix(mglevels[i]->smoothd,&prefix);CHKERRQ(ierr); 1788 ierr = KSPSetOptionsPrefix(subksp,prefix);CHKERRQ(ierr); 1789 ierr = KSPAppendOptionsPrefix(subksp,"up_");CHKERRQ(ierr); 1790 } 1791 PetscFunctionReturn(0); 1792 } 1793 1794 /* No new matrices are created, and the coarse operator matrices are the references to the original ones */ 1795 PetscErrorCode PCGetInterpolations_MG(PC pc,PetscInt *num_levels,Mat *interpolations[]) 1796 { 1797 PC_MG *mg = (PC_MG*)pc->data; 1798 PC_MG_Levels **mglevels = mg->levels; 1799 Mat *mat; 1800 PetscInt l; 1801 PetscErrorCode ierr; 1802 1803 PetscFunctionBegin; 1804 PetscCheckFalse(!mglevels,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling"); 1805 ierr = PetscMalloc1(mg->nlevels,&mat);CHKERRQ(ierr); 1806 for (l=1; l< mg->nlevels; l++) { 1807 mat[l-1] = mglevels[l]->interpolate; 1808 ierr = PetscObjectReference((PetscObject)mat[l-1]);CHKERRQ(ierr); 1809 } 1810 *num_levels = mg->nlevels; 1811 *interpolations = mat; 1812 PetscFunctionReturn(0); 1813 } 1814 1815 /* No new matrices are created, and the coarse operator matrices are the references to the original ones */ 1816 PetscErrorCode PCGetCoarseOperators_MG(PC pc,PetscInt *num_levels,Mat *coarseOperators[]) 1817 { 1818 PC_MG *mg = (PC_MG*)pc->data; 1819 PC_MG_Levels **mglevels = mg->levels; 1820 PetscInt l; 1821 Mat *mat; 1822 PetscErrorCode ierr; 1823 1824 PetscFunctionBegin; 1825 PetscCheckFalse(!mglevels,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling"); 1826 ierr = PetscMalloc1(mg->nlevels,&mat);CHKERRQ(ierr); 1827 for (l=0; l<mg->nlevels-1; l++) { 1828 ierr = KSPGetOperators(mglevels[l]->smoothd,NULL,&(mat[l]));CHKERRQ(ierr); 1829 ierr = PetscObjectReference((PetscObject)mat[l]);CHKERRQ(ierr); 1830 } 1831 *num_levels = mg->nlevels; 1832 *coarseOperators = mat; 1833 PetscFunctionReturn(0); 1834 } 1835 1836 /*@C 1837 PCMGRegisterCoarseSpaceConstructor - Adds a method to the PCMG package for coarse space construction. 1838 1839 Not collective 1840 1841 Input Parameters: 1842 + name - name of the constructor 1843 - function - constructor routine 1844 1845 Notes: 1846 Calling sequence for the routine: 1847 $ my_csp(PC pc, PetscInt l, DM dm, KSP smooth, PetscInt Nc, const Vec initGuess[], Vec **coarseSp) 1848 $ pc - The PC object 1849 $ l - The multigrid level, 0 is the coarse level 1850 $ dm - The DM for this level 1851 $ smooth - The level smoother 1852 $ Nc - The size of the coarse space 1853 $ initGuess - Basis for an initial guess for the space 1854 $ coarseSp - A basis for the computed coarse space 1855 1856 Level: advanced 1857 1858 .seealso: PCMGGetCoarseSpaceConstructor(), PCRegister() 1859 @*/ 1860 PetscErrorCode PCMGRegisterCoarseSpaceConstructor(const char name[], PetscErrorCode (*function)(PC, PetscInt, DM, KSP, PetscInt, const Vec[], Vec **)) 1861 { 1862 PetscErrorCode ierr; 1863 1864 PetscFunctionBegin; 1865 ierr = PCInitializePackage();CHKERRQ(ierr); 1866 ierr = PetscFunctionListAdd(&PCMGCoarseList,name,function);CHKERRQ(ierr); 1867 PetscFunctionReturn(0); 1868 } 1869 1870 /*@C 1871 PCMGGetCoarseSpaceConstructor - Returns the given coarse space construction method. 1872 1873 Not collective 1874 1875 Input Parameter: 1876 . name - name of the constructor 1877 1878 Output Parameter: 1879 . function - constructor routine 1880 1881 Notes: 1882 Calling sequence for the routine: 1883 $ my_csp(PC pc, PetscInt l, DM dm, KSP smooth, PetscInt Nc, const Vec initGuess[], Vec **coarseSp) 1884 $ pc - The PC object 1885 $ l - The multigrid level, 0 is the coarse level 1886 $ dm - The DM for this level 1887 $ smooth - The level smoother 1888 $ Nc - The size of the coarse space 1889 $ initGuess - Basis for an initial guess for the space 1890 $ coarseSp - A basis for the computed coarse space 1891 1892 Level: advanced 1893 1894 .seealso: PCMGRegisterCoarseSpaceConstructor(), PCRegister() 1895 @*/ 1896 PetscErrorCode PCMGGetCoarseSpaceConstructor(const char name[], PetscErrorCode (**function)(PC, PetscInt, DM, KSP, PetscInt, const Vec[], Vec **)) 1897 { 1898 PetscErrorCode ierr; 1899 1900 PetscFunctionBegin; 1901 ierr = PetscFunctionListFind(PCMGCoarseList,name,function);CHKERRQ(ierr); 1902 PetscFunctionReturn(0); 1903 } 1904 1905 /* ----------------------------------------------------------------------------------------*/ 1906 1907 /*MC 1908 PCMG - Use multigrid preconditioning. This preconditioner requires you provide additional 1909 information about the coarser grid matrices and restriction/interpolation operators. 1910 1911 Options Database Keys: 1912 + -pc_mg_levels <nlevels> - number of levels including finest 1913 . -pc_mg_cycle_type <v,w> - provide the cycle desired 1914 . -pc_mg_type <additive,multiplicative,full,kaskade> - multiplicative is the default 1915 . -pc_mg_log - log information about time spent on each level of the solver 1916 . -pc_mg_distinct_smoothup - configure up (after interpolation) and down (before restriction) smoothers separately (with different options prefixes) 1917 . -pc_mg_galerkin <both,pmat,mat,none> - use Galerkin process to compute coarser operators, i.e. Acoarse = R A R' 1918 . -pc_mg_multiplicative_cycles - number of cycles to use as the preconditioner (defaults to 1) 1919 . -pc_mg_dump_matlab - dumps the matrices for each level and the restriction/interpolation matrices 1920 to the Socket viewer for reading from MATLAB. 1921 - -pc_mg_dump_binary - dumps the matrices for each level and the restriction/interpolation matrices 1922 to the binary output file called binaryoutput 1923 1924 Notes: 1925 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 1926 1927 When run with a single level the smoother options are used on that level NOT the coarse grid solver options 1928 1929 When run with KSPRICHARDSON the convergence test changes slightly if monitor is turned on. The iteration count may change slightly. This 1930 is because without monitoring the residual norm is computed WITHIN each multigrid cycle on the finest level after the pre-smoothing 1931 (because the residual has just been computed for the multigrid algorithm and is hence available for free) while with monitoring the 1932 residual is computed at the end of each cycle. 1933 1934 Level: intermediate 1935 1936 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCMGType, PCEXOTIC, PCGAMG, PCML, PCHYPRE 1937 PCMGSetLevels(), PCMGGetLevels(), PCMGSetType(), PCMGSetCycleType(), 1938 PCMGSetDistinctSmoothUp(), PCMGGetCoarseSolve(), PCMGSetResidual(), PCMGSetInterpolation(), 1939 PCMGSetRestriction(), PCMGGetSmoother(), PCMGGetSmootherUp(), PCMGGetSmootherDown(), 1940 PCMGSetCycleTypeOnLevel(), PCMGSetRhs(), PCMGSetX(), PCMGSetR() 1941 M*/ 1942 1943 PETSC_EXTERN PetscErrorCode PCCreate_MG(PC pc) 1944 { 1945 PC_MG *mg; 1946 PetscErrorCode ierr; 1947 1948 PetscFunctionBegin; 1949 ierr = PetscNewLog(pc,&mg);CHKERRQ(ierr); 1950 pc->data = mg; 1951 mg->nlevels = -1; 1952 mg->am = PC_MG_MULTIPLICATIVE; 1953 mg->galerkin = PC_MG_GALERKIN_NONE; 1954 mg->adaptInterpolation = PETSC_FALSE; 1955 mg->Nc = -1; 1956 mg->eigenvalue = -1; 1957 1958 pc->useAmat = PETSC_TRUE; 1959 1960 pc->ops->apply = PCApply_MG; 1961 pc->ops->applytranspose = PCApplyTranspose_MG; 1962 pc->ops->matapply = PCMatApply_MG; 1963 pc->ops->setup = PCSetUp_MG; 1964 pc->ops->reset = PCReset_MG; 1965 pc->ops->destroy = PCDestroy_MG; 1966 pc->ops->setfromoptions = PCSetFromOptions_MG; 1967 pc->ops->view = PCView_MG; 1968 1969 ierr = PetscObjectComposedDataRegister(&mg->eigenvalue);CHKERRQ(ierr); 1970 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCMGSetGalerkin_C",PCMGSetGalerkin_MG);CHKERRQ(ierr); 1971 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCMGGetLevels_C",PCMGGetLevels_MG);CHKERRQ(ierr); 1972 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCMGSetLevels_C",PCMGSetLevels_MG);CHKERRQ(ierr); 1973 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGetInterpolations_C",PCGetInterpolations_MG);CHKERRQ(ierr); 1974 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGetCoarseOperators_C",PCGetCoarseOperators_MG);CHKERRQ(ierr); 1975 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCMGSetAdaptInterpolation_C",PCMGSetAdaptInterpolation_MG);CHKERRQ(ierr); 1976 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCMGGetAdaptInterpolation_C",PCMGGetAdaptInterpolation_MG);CHKERRQ(ierr); 1977 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCMGSetAdaptCR_C",PCMGSetAdaptCR_MG);CHKERRQ(ierr); 1978 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCMGGetAdaptCR_C",PCMGGetAdaptCR_MG);CHKERRQ(ierr); 1979 PetscFunctionReturn(0); 1980 } 1981