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