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, (void **) &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, (void **) &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, (void **) &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. 463 464 Level: intermediate 465 466 Notes: 467 If the number of levels is one then the multigrid uses the -mg_levels prefix 468 for setting the level options rather than the -mg_coarse prefix. 469 470 You can free the information in comms after this routine is called. 471 472 The array of MPI communicators must contain MPI_COMM_NULL for those ranks that at each level 473 are not participating in the coarser solve. For example, with 2 levels and 1 and 2 ranks on 474 the two levels, rank 0 in the original communicator will pass in an array of 2 communicators 475 of size 2 and 1, while rank 1 in the original communicator will pass in array of 2 communicators 476 the first of size 2 and the second of value MPI_COMM_NULL since the rank 1 does not participate 477 in the coarse grid solve. 478 479 Since each coarser level may have a new MPI_Comm with fewer ranks than the previous, one 480 must take special care in providing the restriction and interpolation operation. We recommend 481 providing these as two step operations; first perform a standard restriction or interpolation on 482 the full number of ranks for that level and then use an MPI call to copy the resulting vector 483 array entries (after calls to VecGetArray()) to the smaller or larger number of ranks, not in both 484 cases the MPI calls must be made on the larger of the two communicators. Traditional MPI send and 485 recieves or MPI_AlltoAllv() could be used to do the reshuffling of the vector entries. 486 487 488 .seealso: PCMGSetType(), PCMGGetLevels() 489 @*/ 490 PetscErrorCode PCMGSetLevels(PC pc,PetscInt levels,MPI_Comm *comms) 491 { 492 PetscErrorCode ierr; 493 494 PetscFunctionBegin; 495 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 496 if (comms) PetscValidPointer(comms,3); 497 ierr = PetscTryMethod(pc,"PCMGSetLevels_C",(PC,PetscInt,MPI_Comm*),(pc,levels,comms));CHKERRQ(ierr); 498 PetscFunctionReturn(0); 499 } 500 501 502 PetscErrorCode PCDestroy_MG(PC pc) 503 { 504 PetscErrorCode ierr; 505 PC_MG *mg = (PC_MG*)pc->data; 506 PC_MG_Levels **mglevels = mg->levels; 507 PetscInt i,n; 508 509 PetscFunctionBegin; 510 ierr = PCReset_MG(pc);CHKERRQ(ierr); 511 if (mglevels) { 512 n = mglevels[0]->levels; 513 for (i=0; i<n; i++) { 514 if (mglevels[i]->smoothd != mglevels[i]->smoothu) { 515 ierr = KSPDestroy(&mglevels[i]->smoothd);CHKERRQ(ierr); 516 } 517 ierr = KSPDestroy(&mglevels[i]->smoothu);CHKERRQ(ierr); 518 ierr = KSPDestroy(&mglevels[i]->cr);CHKERRQ(ierr); 519 ierr = PetscFree(mglevels[i]);CHKERRQ(ierr); 520 } 521 ierr = PetscFree(mg->levels);CHKERRQ(ierr); 522 } 523 ierr = PetscFree(pc->data);CHKERRQ(ierr); 524 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGetInterpolations_C",NULL);CHKERRQ(ierr); 525 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGetCoarseOperators_C",NULL);CHKERRQ(ierr); 526 PetscFunctionReturn(0); 527 } 528 529 530 /* 531 PCApply_MG - Runs either an additive, multiplicative, Kaskadic 532 or full cycle of multigrid. 533 534 Note: 535 A simple wrapper which calls PCMGMCycle(),PCMGACycle(), or PCMGFCycle(). 536 */ 537 static PetscErrorCode PCApply_MG_Internal(PC pc,Vec b,Vec x,Mat B,Mat X,PetscBool transpose) 538 { 539 PC_MG *mg = (PC_MG*)pc->data; 540 PC_MG_Levels **mglevels = mg->levels; 541 PetscErrorCode ierr; 542 PC tpc; 543 PetscInt levels = mglevels[0]->levels,i; 544 PetscBool changeu,changed,matapp; 545 546 PetscFunctionBegin; 547 matapp = (PetscBool)(B && X); 548 if (mg->stageApply) {ierr = PetscLogStagePush(mg->stageApply);CHKERRQ(ierr);} 549 /* When the DM is supplying the matrix then it will not exist until here */ 550 for (i=0; i<levels; i++) { 551 if (!mglevels[i]->A) { 552 ierr = KSPGetOperators(mglevels[i]->smoothu,&mglevels[i]->A,NULL);CHKERRQ(ierr); 553 ierr = PetscObjectReference((PetscObject)mglevels[i]->A);CHKERRQ(ierr); 554 } 555 } 556 557 ierr = KSPGetPC(mglevels[levels-1]->smoothd,&tpc);CHKERRQ(ierr); 558 ierr = PCPreSolveChangeRHS(tpc,&changed);CHKERRQ(ierr); 559 ierr = KSPGetPC(mglevels[levels-1]->smoothu,&tpc);CHKERRQ(ierr); 560 ierr = PCPreSolveChangeRHS(tpc,&changeu);CHKERRQ(ierr); 561 if (!changeu && !changed) { 562 if (matapp) { 563 ierr = MatDestroy(&mglevels[levels-1]->B);CHKERRQ(ierr); 564 mglevels[levels-1]->B = B; 565 } else { 566 ierr = VecDestroy(&mglevels[levels-1]->b);CHKERRQ(ierr); 567 mglevels[levels-1]->b = b; 568 } 569 } else { /* if the smoother changes the rhs during PreSolve, we cannot use the input vector */ 570 if (matapp) { 571 if (mglevels[levels-1]->B) { 572 PetscInt N1,N2; 573 PetscBool flg; 574 575 ierr = MatGetSize(mglevels[levels-1]->B,NULL,&N1);CHKERRQ(ierr); 576 ierr = MatGetSize(B,NULL,&N2);CHKERRQ(ierr); 577 ierr = PetscObjectTypeCompare((PetscObject)mglevels[levels-1]->B,((PetscObject)B)->type_name,&flg);CHKERRQ(ierr); 578 if (N1 != N2 || !flg) { 579 ierr = MatDestroy(&mglevels[levels-1]->B);CHKERRQ(ierr); 580 } 581 } 582 if (!mglevels[levels-1]->B) { 583 ierr = MatDuplicate(B,MAT_COPY_VALUES,&mglevels[levels-1]->B);CHKERRQ(ierr); 584 } else { 585 ierr = MatCopy(B,mglevels[levels-1]->B,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 586 } 587 } else { 588 if (!mglevels[levels-1]->b) { 589 Vec *vec; 590 591 ierr = KSPCreateVecs(mglevels[levels-1]->smoothd,1,&vec,0,NULL);CHKERRQ(ierr); 592 mglevels[levels-1]->b = *vec; 593 ierr = PetscFree(vec);CHKERRQ(ierr); 594 } 595 ierr = VecCopy(b,mglevels[levels-1]->b);CHKERRQ(ierr); 596 } 597 } 598 if (matapp) { mglevels[levels-1]->X = X; } 599 else { mglevels[levels-1]->x = x; } 600 601 /* If coarser Xs are present, it means we have already block applied the PC at least once 602 Reset operators if sizes/type do no match */ 603 if (matapp && levels > 1 && mglevels[levels-2]->X) { 604 PetscInt Xc,Bc; 605 PetscBool flg; 606 607 ierr = MatGetSize(mglevels[levels-2]->X,NULL,&Xc);CHKERRQ(ierr); 608 ierr = MatGetSize(mglevels[levels-1]->B,NULL,&Bc);CHKERRQ(ierr); 609 ierr = PetscObjectTypeCompare((PetscObject)mglevels[levels-2]->X,((PetscObject)mglevels[levels-1]->X)->type_name,&flg);CHKERRQ(ierr); 610 if (Xc != Bc || !flg) { 611 ierr = MatDestroy(&mglevels[levels-1]->R);CHKERRQ(ierr); 612 for (i=0;i<levels-1;i++) { 613 ierr = MatDestroy(&mglevels[i]->R);CHKERRQ(ierr); 614 ierr = MatDestroy(&mglevels[i]->B);CHKERRQ(ierr); 615 ierr = MatDestroy(&mglevels[i]->X);CHKERRQ(ierr); 616 } 617 } 618 } 619 620 if (mg->am == PC_MG_MULTIPLICATIVE) { 621 if (matapp) { ierr = MatZeroEntries(X);CHKERRQ(ierr); } 622 else { ierr = VecZeroEntries(x);CHKERRQ(ierr); } 623 for (i=0; i<mg->cyclesperpcapply; i++) { 624 ierr = PCMGMCycle_Private(pc,mglevels+levels-1,transpose,matapp,NULL);CHKERRQ(ierr); 625 } 626 } else if (mg->am == PC_MG_ADDITIVE) { 627 ierr = PCMGACycle_Private(pc,mglevels,transpose,matapp);CHKERRQ(ierr); 628 } else if (mg->am == PC_MG_KASKADE) { 629 ierr = PCMGKCycle_Private(pc,mglevels,transpose,matapp);CHKERRQ(ierr); 630 } else { 631 ierr = PCMGFCycle_Private(pc,mglevels,transpose,matapp);CHKERRQ(ierr); 632 } 633 if (mg->stageApply) {ierr = PetscLogStagePop();CHKERRQ(ierr);} 634 if (!changeu && !changed) { 635 if (matapp) { mglevels[levels-1]->B = NULL; } 636 else { mglevels[levels-1]->b = NULL; } 637 } 638 PetscFunctionReturn(0); 639 } 640 641 static PetscErrorCode PCApply_MG(PC pc,Vec b,Vec x) 642 { 643 PetscErrorCode ierr; 644 645 PetscFunctionBegin; 646 ierr = PCApply_MG_Internal(pc,b,x,NULL,NULL,PETSC_FALSE);CHKERRQ(ierr); 647 PetscFunctionReturn(0); 648 } 649 650 static PetscErrorCode PCApplyTranspose_MG(PC pc,Vec b,Vec x) 651 { 652 PetscErrorCode ierr; 653 654 PetscFunctionBegin; 655 ierr = PCApply_MG_Internal(pc,b,x,NULL,NULL,PETSC_TRUE);CHKERRQ(ierr); 656 PetscFunctionReturn(0); 657 } 658 659 static PetscErrorCode PCMatApply_MG(PC pc,Mat b,Mat x) 660 { 661 PetscErrorCode ierr; 662 663 PetscFunctionBegin; 664 ierr = PCApply_MG_Internal(pc,NULL,NULL,b,x,PETSC_FALSE);CHKERRQ(ierr); 665 PetscFunctionReturn(0); 666 } 667 668 PetscErrorCode PCSetFromOptions_MG(PetscOptionItems *PetscOptionsObject,PC pc) 669 { 670 PetscErrorCode ierr; 671 PetscInt levels,cycles; 672 PetscBool flg, flg2; 673 PC_MG *mg = (PC_MG*)pc->data; 674 PC_MG_Levels **mglevels; 675 PCMGType mgtype; 676 PCMGCycleType mgctype; 677 PCMGGalerkinType gtype; 678 679 PetscFunctionBegin; 680 levels = PetscMax(mg->nlevels,1); 681 ierr = PetscOptionsHead(PetscOptionsObject,"Multigrid options");CHKERRQ(ierr); 682 ierr = PetscOptionsInt("-pc_mg_levels","Number of Levels","PCMGSetLevels",levels,&levels,&flg);CHKERRQ(ierr); 683 if (!flg && !mg->levels && pc->dm) { 684 ierr = DMGetRefineLevel(pc->dm,&levels);CHKERRQ(ierr); 685 levels++; 686 mg->usedmfornumberoflevels = PETSC_TRUE; 687 } 688 ierr = PCMGSetLevels(pc,levels,NULL);CHKERRQ(ierr); 689 mglevels = mg->levels; 690 691 mgctype = (PCMGCycleType) mglevels[0]->cycles; 692 ierr = PetscOptionsEnum("-pc_mg_cycle_type","V cycle or for W-cycle","PCMGSetCycleType",PCMGCycleTypes,(PetscEnum)mgctype,(PetscEnum*)&mgctype,&flg);CHKERRQ(ierr); 693 if (flg) { 694 ierr = PCMGSetCycleType(pc,mgctype);CHKERRQ(ierr); 695 } 696 gtype = mg->galerkin; 697 ierr = PetscOptionsEnum("-pc_mg_galerkin","Use Galerkin process to compute coarser operators","PCMGSetGalerkin",PCMGGalerkinTypes,(PetscEnum)gtype,(PetscEnum*)>ype,&flg);CHKERRQ(ierr); 698 if (flg) { 699 ierr = PCMGSetGalerkin(pc,gtype);CHKERRQ(ierr); 700 } 701 flg2 = PETSC_FALSE; 702 ierr = PetscOptionsBool("-pc_mg_adapt_interp","Adapt interpolation using some coarse space","PCMGSetAdaptInterpolation",PETSC_FALSE,&flg2,&flg);CHKERRQ(ierr); 703 if (flg) {ierr = PCMGSetAdaptInterpolation(pc, flg2);CHKERRQ(ierr);} 704 ierr = PetscOptionsInt("-pc_mg_adapt_interp_n","Size of the coarse space for adaptive interpolation","PCMGSetCoarseSpace",mg->Nc,&mg->Nc,&flg);CHKERRQ(ierr); 705 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); 706 ierr = PetscOptionsBool("-pc_mg_mesp_monitor","Monitor the multilevel eigensolver","PCMGSetAdaptInterpolation",PETSC_FALSE,&mg->mespMonitor,&flg);CHKERRQ(ierr); 707 flg2 = PETSC_FALSE; 708 ierr = PetscOptionsBool("-pc_mg_adapt_cr","Monitor coarse space quality using Compatible Relaxation (CR)","PCMGSetAdaptCR",PETSC_FALSE,&flg2,&flg);CHKERRQ(ierr); 709 if (flg) {ierr = PCMGSetAdaptCR(pc, flg2);CHKERRQ(ierr);} 710 flg = PETSC_FALSE; 711 ierr = PetscOptionsBool("-pc_mg_distinct_smoothup","Create separate smoothup KSP and append the prefix _up","PCMGSetDistinctSmoothUp",PETSC_FALSE,&flg,NULL);CHKERRQ(ierr); 712 if (flg) { 713 ierr = PCMGSetDistinctSmoothUp(pc);CHKERRQ(ierr); 714 } 715 mgtype = mg->am; 716 ierr = PetscOptionsEnum("-pc_mg_type","Multigrid type","PCMGSetType",PCMGTypes,(PetscEnum)mgtype,(PetscEnum*)&mgtype,&flg);CHKERRQ(ierr); 717 if (flg) { 718 ierr = PCMGSetType(pc,mgtype);CHKERRQ(ierr); 719 } 720 if (mg->am == PC_MG_MULTIPLICATIVE) { 721 ierr = PetscOptionsInt("-pc_mg_multiplicative_cycles","Number of cycles for each preconditioner step","PCMGMultiplicativeSetCycles",mg->cyclesperpcapply,&cycles,&flg);CHKERRQ(ierr); 722 if (flg) { 723 ierr = PCMGMultiplicativeSetCycles(pc,cycles);CHKERRQ(ierr); 724 } 725 } 726 flg = PETSC_FALSE; 727 ierr = PetscOptionsBool("-pc_mg_log","Log times for each multigrid level","None",flg,&flg,NULL);CHKERRQ(ierr); 728 if (flg) { 729 PetscInt i; 730 char eventname[128]; 731 732 levels = mglevels[0]->levels; 733 for (i=0; i<levels; i++) { 734 sprintf(eventname,"MGSetup Level %d",(int)i); 735 ierr = PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventsmoothsetup);CHKERRQ(ierr); 736 sprintf(eventname,"MGSmooth Level %d",(int)i); 737 ierr = PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventsmoothsolve);CHKERRQ(ierr); 738 if (i) { 739 sprintf(eventname,"MGResid Level %d",(int)i); 740 ierr = PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventresidual);CHKERRQ(ierr); 741 sprintf(eventname,"MGInterp Level %d",(int)i); 742 ierr = PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventinterprestrict);CHKERRQ(ierr); 743 } 744 } 745 746 #if defined(PETSC_USE_LOG) 747 { 748 const char *sname = "MG Apply"; 749 PetscStageLog stageLog; 750 PetscInt st; 751 752 ierr = PetscLogGetStageLog(&stageLog);CHKERRQ(ierr); 753 for (st = 0; st < stageLog->numStages; ++st) { 754 PetscBool same; 755 756 ierr = PetscStrcmp(stageLog->stageInfo[st].name, sname, &same);CHKERRQ(ierr); 757 if (same) mg->stageApply = st; 758 } 759 if (!mg->stageApply) { 760 ierr = PetscLogStageRegister(sname, &mg->stageApply);CHKERRQ(ierr); 761 } 762 } 763 #endif 764 } 765 ierr = PetscOptionsTail();CHKERRQ(ierr); 766 /* Check option consistency */ 767 ierr = PCMGGetGalerkin(pc, >ype);CHKERRQ(ierr); 768 ierr = PCMGGetAdaptInterpolation(pc, &flg);CHKERRQ(ierr); 769 if (flg && (gtype >= PC_MG_GALERKIN_NONE)) SETERRQ(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_INCOMP, "Must use Galerkin coarse operators when adapting the interpolator"); 770 PetscFunctionReturn(0); 771 } 772 773 const char *const PCMGTypes[] = {"MULTIPLICATIVE","ADDITIVE","FULL","KASKADE","PCMGType","PC_MG",NULL}; 774 const char *const PCMGCycleTypes[] = {"invalid","v","w","PCMGCycleType","PC_MG_CYCLE",NULL}; 775 const char *const PCMGGalerkinTypes[] = {"both","pmat","mat","none","external","PCMGGalerkinType","PC_MG_GALERKIN",NULL}; 776 const char *const PCMGCoarseSpaceTypes[] = {"polynomial","harmonic","eigenvector","generalized_eigenvector","PCMGCoarseSpaceType","PCMG_POLYNOMIAL",NULL}; 777 778 #include <petscdraw.h> 779 PetscErrorCode PCView_MG(PC pc,PetscViewer viewer) 780 { 781 PC_MG *mg = (PC_MG*)pc->data; 782 PC_MG_Levels **mglevels = mg->levels; 783 PetscErrorCode ierr; 784 PetscInt levels = mglevels ? mglevels[0]->levels : 0,i; 785 PetscBool iascii,isbinary,isdraw; 786 787 PetscFunctionBegin; 788 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 789 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 790 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 791 if (iascii) { 792 const char *cyclename = levels ? (mglevels[0]->cycles == PC_MG_CYCLE_V ? "v" : "w") : "unknown"; 793 ierr = PetscViewerASCIIPrintf(viewer," type is %s, levels=%D cycles=%s\n", PCMGTypes[mg->am],levels,cyclename);CHKERRQ(ierr); 794 if (mg->am == PC_MG_MULTIPLICATIVE) { 795 ierr = PetscViewerASCIIPrintf(viewer," Cycles per PCApply=%d\n",mg->cyclesperpcapply);CHKERRQ(ierr); 796 } 797 if (mg->galerkin == PC_MG_GALERKIN_BOTH) { 798 ierr = PetscViewerASCIIPrintf(viewer," Using Galerkin computed coarse grid matrices\n");CHKERRQ(ierr); 799 } else if (mg->galerkin == PC_MG_GALERKIN_PMAT) { 800 ierr = PetscViewerASCIIPrintf(viewer," Using Galerkin computed coarse grid matrices for pmat\n");CHKERRQ(ierr); 801 } else if (mg->galerkin == PC_MG_GALERKIN_MAT) { 802 ierr = PetscViewerASCIIPrintf(viewer," Using Galerkin computed coarse grid matrices for mat\n");CHKERRQ(ierr); 803 } else if (mg->galerkin == PC_MG_GALERKIN_EXTERNAL) { 804 ierr = PetscViewerASCIIPrintf(viewer," Using externally compute Galerkin coarse grid matrices\n");CHKERRQ(ierr); 805 } else { 806 ierr = PetscViewerASCIIPrintf(viewer," Not using Galerkin computed coarse grid matrices\n");CHKERRQ(ierr); 807 } 808 if (mg->view){ 809 ierr = (*mg->view)(pc,viewer);CHKERRQ(ierr); 810 } 811 for (i=0; i<levels; i++) { 812 if (!i) { 813 ierr = PetscViewerASCIIPrintf(viewer,"Coarse grid solver -- level -------------------------------\n",i);CHKERRQ(ierr); 814 } else { 815 ierr = PetscViewerASCIIPrintf(viewer,"Down solver (pre-smoother) on level %D -------------------------------\n",i);CHKERRQ(ierr); 816 } 817 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 818 ierr = KSPView(mglevels[i]->smoothd,viewer);CHKERRQ(ierr); 819 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 820 if (i && mglevels[i]->smoothd == mglevels[i]->smoothu) { 821 ierr = PetscViewerASCIIPrintf(viewer,"Up solver (post-smoother) same as down solver (pre-smoother)\n");CHKERRQ(ierr); 822 } else if (i) { 823 ierr = PetscViewerASCIIPrintf(viewer,"Up solver (post-smoother) on level %D -------------------------------\n",i);CHKERRQ(ierr); 824 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 825 ierr = KSPView(mglevels[i]->smoothu,viewer);CHKERRQ(ierr); 826 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 827 } 828 if (i && mglevels[i]->cr) { 829 ierr = PetscViewerASCIIPrintf(viewer,"CR solver on level %D -------------------------------\n",i);CHKERRQ(ierr); 830 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 831 ierr = KSPView(mglevels[i]->cr,viewer);CHKERRQ(ierr); 832 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 833 } 834 } 835 } else if (isbinary) { 836 for (i=levels-1; i>=0; i--) { 837 ierr = KSPView(mglevels[i]->smoothd,viewer);CHKERRQ(ierr); 838 if (i && mglevels[i]->smoothd != mglevels[i]->smoothu) { 839 ierr = KSPView(mglevels[i]->smoothu,viewer);CHKERRQ(ierr); 840 } 841 } 842 } else if (isdraw) { 843 PetscDraw draw; 844 PetscReal x,w,y,bottom,th; 845 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 846 ierr = PetscDrawGetCurrentPoint(draw,&x,&y);CHKERRQ(ierr); 847 ierr = PetscDrawStringGetSize(draw,NULL,&th);CHKERRQ(ierr); 848 bottom = y - th; 849 for (i=levels-1; i>=0; i--) { 850 if (!mglevels[i]->smoothu || (mglevels[i]->smoothu == mglevels[i]->smoothd)) { 851 ierr = PetscDrawPushCurrentPoint(draw,x,bottom);CHKERRQ(ierr); 852 ierr = KSPView(mglevels[i]->smoothd,viewer);CHKERRQ(ierr); 853 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 854 } else { 855 w = 0.5*PetscMin(1.0-x,x); 856 ierr = PetscDrawPushCurrentPoint(draw,x+w,bottom);CHKERRQ(ierr); 857 ierr = KSPView(mglevels[i]->smoothd,viewer);CHKERRQ(ierr); 858 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 859 ierr = PetscDrawPushCurrentPoint(draw,x-w,bottom);CHKERRQ(ierr); 860 ierr = KSPView(mglevels[i]->smoothu,viewer);CHKERRQ(ierr); 861 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 862 } 863 ierr = PetscDrawGetBoundingBox(draw,NULL,&bottom,NULL,NULL);CHKERRQ(ierr); 864 bottom -= th; 865 } 866 } 867 PetscFunctionReturn(0); 868 } 869 870 #include <petsc/private/kspimpl.h> 871 872 /* 873 Calls setup for the KSP on each level 874 */ 875 PetscErrorCode PCSetUp_MG(PC pc) 876 { 877 PC_MG *mg = (PC_MG*)pc->data; 878 PC_MG_Levels **mglevels = mg->levels; 879 PetscErrorCode ierr; 880 PetscInt i,n; 881 PC cpc; 882 PetscBool dump = PETSC_FALSE,opsset,use_amat,missinginterpolate = PETSC_FALSE; 883 Mat dA,dB; 884 Vec tvec; 885 DM *dms; 886 PetscViewer viewer = NULL; 887 PetscBool dAeqdB = PETSC_FALSE, needRestricts = PETSC_FALSE, doCR = PETSC_FALSE; 888 889 PetscFunctionBegin; 890 if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels with PCMGSetLevels() before setting up"); 891 n = mglevels[0]->levels; 892 /* FIX: Move this to PCSetFromOptions_MG? */ 893 if (mg->usedmfornumberoflevels) { 894 PetscInt levels; 895 ierr = DMGetRefineLevel(pc->dm,&levels);CHKERRQ(ierr); 896 levels++; 897 if (levels > n) { /* the problem is now being solved on a finer grid */ 898 ierr = PCMGSetLevels(pc,levels,NULL);CHKERRQ(ierr); 899 n = levels; 900 ierr = PCSetFromOptions(pc);CHKERRQ(ierr); /* it is bad to call this here, but otherwise will never be called for the new hierarchy */ 901 mglevels = mg->levels; 902 } 903 } 904 ierr = KSPGetPC(mglevels[0]->smoothd,&cpc);CHKERRQ(ierr); 905 906 907 /* If user did not provide fine grid operators OR operator was not updated since last global KSPSetOperators() */ 908 /* so use those from global PC */ 909 /* Is this what we always want? What if user wants to keep old one? */ 910 ierr = KSPGetOperatorsSet(mglevels[n-1]->smoothd,NULL,&opsset);CHKERRQ(ierr); 911 if (opsset) { 912 Mat mmat; 913 ierr = KSPGetOperators(mglevels[n-1]->smoothd,NULL,&mmat);CHKERRQ(ierr); 914 if (mmat == pc->pmat) opsset = PETSC_FALSE; 915 } 916 917 /* Create CR solvers */ 918 ierr = PCMGGetAdaptCR(pc, &doCR);CHKERRQ(ierr); 919 if (doCR) { 920 const char *prefix; 921 922 ierr = PCGetOptionsPrefix(pc, &prefix);CHKERRQ(ierr); 923 for (i = 1; i < n; ++i) { 924 PC ipc, cr; 925 char crprefix[128]; 926 927 ierr = KSPCreate(PetscObjectComm((PetscObject) pc), &mglevels[i]->cr);CHKERRQ(ierr); 928 ierr = KSPSetErrorIfNotConverged(mglevels[i]->cr, PETSC_FALSE);CHKERRQ(ierr); 929 ierr = PetscObjectIncrementTabLevel((PetscObject) mglevels[i]->cr, (PetscObject) pc, n-i);CHKERRQ(ierr); 930 ierr = KSPSetOptionsPrefix(mglevels[i]->cr, prefix);CHKERRQ(ierr); 931 ierr = PetscObjectComposedDataSetInt((PetscObject) mglevels[i]->cr, PetscMGLevelId, mglevels[i]->level);CHKERRQ(ierr); 932 ierr = KSPSetType(mglevels[i]->cr, KSPCHEBYSHEV);CHKERRQ(ierr); 933 ierr = KSPSetConvergenceTest(mglevels[i]->cr, KSPConvergedSkip, NULL, NULL);CHKERRQ(ierr); 934 ierr = KSPSetNormType(mglevels[i]->cr, KSP_NORM_PRECONDITIONED);CHKERRQ(ierr); 935 ierr = KSPGetPC(mglevels[i]->cr, &ipc);CHKERRQ(ierr); 936 937 ierr = PCSetType(ipc, PCCOMPOSITE);CHKERRQ(ierr); 938 ierr = PCCompositeSetType(ipc, PC_COMPOSITE_MULTIPLICATIVE);CHKERRQ(ierr); 939 ierr = PCCompositeAddPCType(ipc, PCSOR);CHKERRQ(ierr); 940 ierr = CreateCR_Private(pc, i, &cr);CHKERRQ(ierr); 941 ierr = PCCompositeAddPC(ipc, cr);CHKERRQ(ierr); 942 ierr = PCDestroy(&cr);CHKERRQ(ierr); 943 944 ierr = KSPSetTolerances(mglevels[i]->cr, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, mg->default_smoothd);CHKERRQ(ierr); 945 ierr = KSPSetInitialGuessNonzero(mglevels[i]->cr, PETSC_TRUE);CHKERRQ(ierr); 946 ierr = PetscSNPrintf(crprefix, 128, "mg_levels_%d_cr_", (int) i);CHKERRQ(ierr); 947 ierr = KSPAppendOptionsPrefix(mglevels[i]->cr, crprefix);CHKERRQ(ierr); 948 ierr = PetscLogObjectParent((PetscObject) pc, (PetscObject) mglevels[i]->cr);CHKERRQ(ierr); 949 } 950 } 951 952 if (!opsset) { 953 ierr = PCGetUseAmat(pc,&use_amat);CHKERRQ(ierr); 954 if (use_amat) { 955 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); 956 ierr = KSPSetOperators(mglevels[n-1]->smoothd,pc->mat,pc->pmat);CHKERRQ(ierr); 957 } else { 958 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); 959 ierr = KSPSetOperators(mglevels[n-1]->smoothd,pc->pmat,pc->pmat);CHKERRQ(ierr); 960 } 961 } 962 963 for (i=n-1; i>0; i--) { 964 if (!(mglevels[i]->interpolate || mglevels[i]->restrct)) { 965 missinginterpolate = PETSC_TRUE; 966 continue; 967 } 968 } 969 970 ierr = KSPGetOperators(mglevels[n-1]->smoothd,&dA,&dB);CHKERRQ(ierr); 971 if (dA == dB) dAeqdB = PETSC_TRUE; 972 if ((mg->galerkin == PC_MG_GALERKIN_NONE) || (((mg->galerkin == PC_MG_GALERKIN_PMAT) || (mg->galerkin == PC_MG_GALERKIN_MAT)) && !dAeqdB)) { 973 needRestricts = PETSC_TRUE; /* user must compute either mat, pmat, or both so must restrict x to coarser levels */ 974 } 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 1105 /* Adapt interpolation matrices */ 1106 if (mg->adaptInterpolation) { 1107 mg->Nc = mg->Nc < 0 ? 6 : mg->Nc; /* Default to 6 modes */ 1108 for (i = 0; i < n; ++i) { 1109 ierr = PCMGComputeCoarseSpace_Internal(pc, i, mg->coarseSpaceType, mg->Nc, !i ? NULL : mglevels[i-1]->coarseSpace, &mglevels[i]->coarseSpace);CHKERRQ(ierr); 1110 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);} 1111 } 1112 for (i = n-2; i > -1; --i) { 1113 ierr = PCMGRecomputeLevelOperators_Internal(pc, i);CHKERRQ(ierr); 1114 } 1115 } 1116 1117 if (needRestricts && pc->dm) { 1118 for (i=n-2; i>=0; i--) { 1119 DM dmfine,dmcoarse; 1120 Mat Restrict,Inject; 1121 Vec rscale; 1122 ierr = KSPGetDM(mglevels[i+1]->smoothd,&dmfine);CHKERRQ(ierr); 1123 ierr = KSPGetDM(mglevels[i]->smoothd,&dmcoarse);CHKERRQ(ierr); 1124 ierr = PCMGGetRestriction(pc,i+1,&Restrict);CHKERRQ(ierr); 1125 ierr = PCMGGetRScale(pc,i+1,&rscale);CHKERRQ(ierr); 1126 ierr = PCMGGetInjection(pc,i+1,&Inject);CHKERRQ(ierr); 1127 ierr = DMRestrict(dmfine,Restrict,rscale,Inject,dmcoarse);CHKERRQ(ierr); 1128 } 1129 } 1130 1131 if (!pc->setupcalled) { 1132 for (i=0; i<n; i++) { 1133 ierr = KSPSetFromOptions(mglevels[i]->smoothd);CHKERRQ(ierr); 1134 } 1135 for (i=1; i<n; i++) { 1136 if (mglevels[i]->smoothu && (mglevels[i]->smoothu != mglevels[i]->smoothd)) { 1137 ierr = KSPSetFromOptions(mglevels[i]->smoothu);CHKERRQ(ierr); 1138 } 1139 if (mglevels[i]->cr) { 1140 ierr = KSPSetFromOptions(mglevels[i]->cr);CHKERRQ(ierr); 1141 } 1142 } 1143 /* insure that if either interpolation or restriction is set the other other one is set */ 1144 for (i=1; i<n; i++) { 1145 ierr = PCMGGetInterpolation(pc,i,NULL);CHKERRQ(ierr); 1146 ierr = PCMGGetRestriction(pc,i,NULL);CHKERRQ(ierr); 1147 } 1148 for (i=0; i<n-1; i++) { 1149 if (!mglevels[i]->b) { 1150 Vec *vec; 1151 ierr = KSPCreateVecs(mglevels[i]->smoothd,1,&vec,0,NULL);CHKERRQ(ierr); 1152 ierr = PCMGSetRhs(pc,i,*vec);CHKERRQ(ierr); 1153 ierr = VecDestroy(vec);CHKERRQ(ierr); 1154 ierr = PetscFree(vec);CHKERRQ(ierr); 1155 } 1156 if (!mglevels[i]->r && i) { 1157 ierr = VecDuplicate(mglevels[i]->b,&tvec);CHKERRQ(ierr); 1158 ierr = PCMGSetR(pc,i,tvec);CHKERRQ(ierr); 1159 ierr = VecDestroy(&tvec);CHKERRQ(ierr); 1160 } 1161 if (!mglevels[i]->x) { 1162 ierr = VecDuplicate(mglevels[i]->b,&tvec);CHKERRQ(ierr); 1163 ierr = PCMGSetX(pc,i,tvec);CHKERRQ(ierr); 1164 ierr = VecDestroy(&tvec);CHKERRQ(ierr); 1165 } 1166 if (doCR) { 1167 ierr = VecDuplicate(mglevels[i]->b,&mglevels[i]->crx);CHKERRQ(ierr); 1168 ierr = VecDuplicate(mglevels[i]->b,&mglevels[i]->crb);CHKERRQ(ierr); 1169 ierr = VecZeroEntries(mglevels[i]->crb);CHKERRQ(ierr); 1170 } 1171 } 1172 if (n != 1 && !mglevels[n-1]->r) { 1173 /* PCMGSetR() on the finest level if user did not supply it */ 1174 Vec *vec; 1175 ierr = KSPCreateVecs(mglevels[n-1]->smoothd,1,&vec,0,NULL);CHKERRQ(ierr); 1176 ierr = PCMGSetR(pc,n-1,*vec);CHKERRQ(ierr); 1177 ierr = VecDestroy(vec);CHKERRQ(ierr); 1178 ierr = PetscFree(vec);CHKERRQ(ierr); 1179 } 1180 if (doCR) { 1181 ierr = VecDuplicate(mglevels[n-1]->r, &mglevels[n-1]->crx);CHKERRQ(ierr); 1182 ierr = VecDuplicate(mglevels[n-1]->r, &mglevels[n-1]->crb);CHKERRQ(ierr); 1183 ierr = VecZeroEntries(mglevels[n-1]->crb);CHKERRQ(ierr); 1184 } 1185 } 1186 1187 if (pc->dm) { 1188 /* need to tell all the coarser levels to rebuild the matrix using the DM for that level */ 1189 for (i=0; i<n-1; i++) { 1190 if (mglevels[i]->smoothd->setupstage != KSP_SETUP_NEW) mglevels[i]->smoothd->setupstage = KSP_SETUP_NEWMATRIX; 1191 } 1192 } 1193 1194 for (i=1; i<n; i++) { 1195 if (mglevels[i]->smoothu == mglevels[i]->smoothd || mg->am == PC_MG_FULL || mg->am == PC_MG_KASKADE || mg->cyclesperpcapply > 1){ 1196 /* if doing only down then initial guess is zero */ 1197 ierr = KSPSetInitialGuessNonzero(mglevels[i]->smoothd,PETSC_TRUE);CHKERRQ(ierr); 1198 } 1199 if (mglevels[i]->cr) {ierr = KSPSetInitialGuessNonzero(mglevels[i]->cr,PETSC_TRUE);CHKERRQ(ierr);} 1200 if (mglevels[i]->eventsmoothsetup) {ierr = PetscLogEventBegin(mglevels[i]->eventsmoothsetup,0,0,0,0);CHKERRQ(ierr);} 1201 ierr = KSPSetUp(mglevels[i]->smoothd);CHKERRQ(ierr); 1202 if (mglevels[i]->smoothd->reason == KSP_DIVERGED_PC_FAILED) { 1203 pc->failedreason = PC_SUBPC_ERROR; 1204 } 1205 if (mglevels[i]->eventsmoothsetup) {ierr = PetscLogEventEnd(mglevels[i]->eventsmoothsetup,0,0,0,0);CHKERRQ(ierr);} 1206 if (!mglevels[i]->residual) { 1207 Mat mat; 1208 ierr = KSPGetOperators(mglevels[i]->smoothd,&mat,NULL);CHKERRQ(ierr); 1209 ierr = PCMGSetResidual(pc,i,PCMGResidualDefault,mat);CHKERRQ(ierr); 1210 } 1211 if (!mglevels[i]->residualtranspose) { 1212 Mat mat; 1213 ierr = KSPGetOperators(mglevels[i]->smoothd,&mat,NULL);CHKERRQ(ierr); 1214 ierr = PCMGSetResidualTranspose(pc,i,PCMGResidualTransposeDefault,mat);CHKERRQ(ierr); 1215 } 1216 } 1217 for (i=1; i<n; i++) { 1218 if (mglevels[i]->smoothu && mglevels[i]->smoothu != mglevels[i]->smoothd) { 1219 Mat downmat,downpmat; 1220 1221 /* check if operators have been set for up, if not use down operators to set them */ 1222 ierr = KSPGetOperatorsSet(mglevels[i]->smoothu,&opsset,NULL);CHKERRQ(ierr); 1223 if (!opsset) { 1224 ierr = KSPGetOperators(mglevels[i]->smoothd,&downmat,&downpmat);CHKERRQ(ierr); 1225 ierr = KSPSetOperators(mglevels[i]->smoothu,downmat,downpmat);CHKERRQ(ierr); 1226 } 1227 1228 ierr = KSPSetInitialGuessNonzero(mglevels[i]->smoothu,PETSC_TRUE);CHKERRQ(ierr); 1229 if (mglevels[i]->eventsmoothsetup) {ierr = PetscLogEventBegin(mglevels[i]->eventsmoothsetup,0,0,0,0);CHKERRQ(ierr);} 1230 ierr = KSPSetUp(mglevels[i]->smoothu);CHKERRQ(ierr); 1231 if (mglevels[i]->smoothu->reason == KSP_DIVERGED_PC_FAILED) { 1232 pc->failedreason = PC_SUBPC_ERROR; 1233 } 1234 if (mglevels[i]->eventsmoothsetup) {ierr = PetscLogEventEnd(mglevels[i]->eventsmoothsetup,0,0,0,0);CHKERRQ(ierr);} 1235 } 1236 if (mglevels[i]->cr) { 1237 Mat downmat,downpmat; 1238 1239 /* check if operators have been set for up, if not use down operators to set them */ 1240 ierr = KSPGetOperatorsSet(mglevels[i]->cr,&opsset,NULL);CHKERRQ(ierr); 1241 if (!opsset) { 1242 ierr = KSPGetOperators(mglevels[i]->smoothd,&downmat,&downpmat);CHKERRQ(ierr); 1243 ierr = KSPSetOperators(mglevels[i]->cr,downmat,downpmat);CHKERRQ(ierr); 1244 } 1245 1246 ierr = KSPSetInitialGuessNonzero(mglevels[i]->cr,PETSC_TRUE);CHKERRQ(ierr); 1247 if (mglevels[i]->eventsmoothsetup) {ierr = PetscLogEventBegin(mglevels[i]->eventsmoothsetup,0,0,0,0);CHKERRQ(ierr);} 1248 ierr = KSPSetUp(mglevels[i]->cr);CHKERRQ(ierr); 1249 if (mglevels[i]->cr->reason == KSP_DIVERGED_PC_FAILED) { 1250 pc->failedreason = PC_SUBPC_ERROR; 1251 } 1252 if (mglevels[i]->eventsmoothsetup) {ierr = PetscLogEventEnd(mglevels[i]->eventsmoothsetup,0,0,0,0);CHKERRQ(ierr);} 1253 } 1254 } 1255 1256 if (mglevels[0]->eventsmoothsetup) {ierr = PetscLogEventBegin(mglevels[0]->eventsmoothsetup,0,0,0,0);CHKERRQ(ierr);} 1257 ierr = KSPSetUp(mglevels[0]->smoothd);CHKERRQ(ierr); 1258 if (mglevels[0]->smoothd->reason == KSP_DIVERGED_PC_FAILED) { 1259 pc->failedreason = PC_SUBPC_ERROR; 1260 } 1261 if (mglevels[0]->eventsmoothsetup) {ierr = PetscLogEventEnd(mglevels[0]->eventsmoothsetup,0,0,0,0);CHKERRQ(ierr);} 1262 1263 /* 1264 Dump the interpolation/restriction matrices plus the 1265 Jacobian/stiffness on each level. This allows MATLAB users to 1266 easily check if the Galerkin condition A_c = R A_f R^T is satisfied. 1267 1268 Only support one or the other at the same time. 1269 */ 1270 #if defined(PETSC_USE_SOCKET_VIEWER) 1271 ierr = PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_mg_dump_matlab",&dump,NULL);CHKERRQ(ierr); 1272 if (dump) viewer = PETSC_VIEWER_SOCKET_(PetscObjectComm((PetscObject)pc)); 1273 dump = PETSC_FALSE; 1274 #endif 1275 ierr = PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_mg_dump_binary",&dump,NULL);CHKERRQ(ierr); 1276 if (dump) viewer = PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)pc)); 1277 1278 if (viewer) { 1279 for (i=1; i<n; i++) { 1280 ierr = MatView(mglevels[i]->restrct,viewer);CHKERRQ(ierr); 1281 } 1282 for (i=0; i<n; i++) { 1283 ierr = KSPGetPC(mglevels[i]->smoothd,&pc);CHKERRQ(ierr); 1284 ierr = MatView(pc->mat,viewer);CHKERRQ(ierr); 1285 } 1286 } 1287 PetscFunctionReturn(0); 1288 } 1289 1290 /* -------------------------------------------------------------------------------------*/ 1291 1292 PetscErrorCode PCMGGetLevels_MG(PC pc, PetscInt *levels) 1293 { 1294 PC_MG *mg = (PC_MG *) pc->data; 1295 1296 PetscFunctionBegin; 1297 *levels = mg->nlevels; 1298 PetscFunctionReturn(0); 1299 } 1300 1301 /*@ 1302 PCMGGetLevels - Gets the number of levels to use with MG. 1303 1304 Not Collective 1305 1306 Input Parameter: 1307 . pc - the preconditioner context 1308 1309 Output parameter: 1310 . levels - the number of levels 1311 1312 Level: advanced 1313 1314 .seealso: PCMGSetLevels() 1315 @*/ 1316 PetscErrorCode PCMGGetLevels(PC pc,PetscInt *levels) 1317 { 1318 PetscErrorCode ierr; 1319 1320 PetscFunctionBegin; 1321 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1322 PetscValidIntPointer(levels,2); 1323 *levels = 0; 1324 ierr = PetscTryMethod(pc,"PCMGGetLevels_C",(PC,PetscInt*),(pc,levels));CHKERRQ(ierr); 1325 PetscFunctionReturn(0); 1326 } 1327 1328 /*@ 1329 PCMGSetType - Determines the form of multigrid to use: 1330 multiplicative, additive, full, or the Kaskade algorithm. 1331 1332 Logically Collective on PC 1333 1334 Input Parameters: 1335 + pc - the preconditioner context 1336 - form - multigrid form, one of PC_MG_MULTIPLICATIVE, PC_MG_ADDITIVE, 1337 PC_MG_FULL, PC_MG_KASKADE 1338 1339 Options Database Key: 1340 . -pc_mg_type <form> - Sets <form>, one of multiplicative, 1341 additive, full, kaskade 1342 1343 Level: advanced 1344 1345 .seealso: PCMGSetLevels() 1346 @*/ 1347 PetscErrorCode PCMGSetType(PC pc,PCMGType form) 1348 { 1349 PC_MG *mg = (PC_MG*)pc->data; 1350 1351 PetscFunctionBegin; 1352 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1353 PetscValidLogicalCollectiveEnum(pc,form,2); 1354 mg->am = form; 1355 if (form == PC_MG_MULTIPLICATIVE) pc->ops->applyrichardson = PCApplyRichardson_MG; 1356 else pc->ops->applyrichardson = NULL; 1357 PetscFunctionReturn(0); 1358 } 1359 1360 /*@ 1361 PCMGGetType - Determines the form of multigrid to use: 1362 multiplicative, additive, full, or the Kaskade algorithm. 1363 1364 Logically Collective on PC 1365 1366 Input Parameter: 1367 . pc - the preconditioner context 1368 1369 Output Parameter: 1370 . type - one of PC_MG_MULTIPLICATIVE, PC_MG_ADDITIVE,PC_MG_FULL, PC_MG_KASKADE 1371 1372 1373 Level: advanced 1374 1375 .seealso: PCMGSetLevels() 1376 @*/ 1377 PetscErrorCode PCMGGetType(PC pc,PCMGType *type) 1378 { 1379 PC_MG *mg = (PC_MG*)pc->data; 1380 1381 PetscFunctionBegin; 1382 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1383 *type = mg->am; 1384 PetscFunctionReturn(0); 1385 } 1386 1387 /*@ 1388 PCMGSetCycleType - Sets the type cycles to use. Use PCMGSetCycleTypeOnLevel() for more 1389 complicated cycling. 1390 1391 Logically Collective on PC 1392 1393 Input Parameters: 1394 + pc - the multigrid context 1395 - n - either PC_MG_CYCLE_V or PC_MG_CYCLE_W 1396 1397 Options Database Key: 1398 . -pc_mg_cycle_type <v,w> - provide the cycle desired 1399 1400 Level: advanced 1401 1402 .seealso: PCMGSetCycleTypeOnLevel() 1403 @*/ 1404 PetscErrorCode PCMGSetCycleType(PC pc,PCMGCycleType n) 1405 { 1406 PC_MG *mg = (PC_MG*)pc->data; 1407 PC_MG_Levels **mglevels = mg->levels; 1408 PetscInt i,levels; 1409 1410 PetscFunctionBegin; 1411 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1412 PetscValidLogicalCollectiveEnum(pc,n,2); 1413 if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Must set MG levels with PCMGSetLevels() before calling"); 1414 levels = mglevels[0]->levels; 1415 for (i=0; i<levels; i++) mglevels[i]->cycles = n; 1416 PetscFunctionReturn(0); 1417 } 1418 1419 /*@ 1420 PCMGMultiplicativeSetCycles - Sets the number of cycles to use for each preconditioner step 1421 of multigrid when PCMGType of PC_MG_MULTIPLICATIVE is used 1422 1423 Logically Collective on PC 1424 1425 Input Parameters: 1426 + pc - the multigrid context 1427 - n - number of cycles (default is 1) 1428 1429 Options Database Key: 1430 . -pc_mg_multiplicative_cycles n 1431 1432 Level: advanced 1433 1434 Notes: 1435 This is not associated with setting a v or w cycle, that is set with PCMGSetCycleType() 1436 1437 .seealso: PCMGSetCycleTypeOnLevel(), PCMGSetCycleType() 1438 @*/ 1439 PetscErrorCode PCMGMultiplicativeSetCycles(PC pc,PetscInt n) 1440 { 1441 PC_MG *mg = (PC_MG*)pc->data; 1442 1443 PetscFunctionBegin; 1444 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1445 PetscValidLogicalCollectiveInt(pc,n,2); 1446 mg->cyclesperpcapply = n; 1447 PetscFunctionReturn(0); 1448 } 1449 1450 PetscErrorCode PCMGSetGalerkin_MG(PC pc,PCMGGalerkinType use) 1451 { 1452 PC_MG *mg = (PC_MG*)pc->data; 1453 1454 PetscFunctionBegin; 1455 mg->galerkin = use; 1456 PetscFunctionReturn(0); 1457 } 1458 1459 /*@ 1460 PCMGSetGalerkin - Causes the coarser grid matrices to be computed from the 1461 finest grid via the Galerkin process: A_i-1 = r_i * A_i * p_i 1462 1463 Logically Collective on PC 1464 1465 Input Parameters: 1466 + pc - the multigrid context 1467 - use - one of PC_MG_GALERKIN_BOTH,PC_MG_GALERKIN_PMAT,PC_MG_GALERKIN_MAT, or PC_MG_GALERKIN_NONE 1468 1469 Options Database Key: 1470 . -pc_mg_galerkin <both,pmat,mat,none> 1471 1472 Level: intermediate 1473 1474 Notes: 1475 Some codes that use PCMG such as PCGAMG use Galerkin internally while constructing the hierarchy and thus do not 1476 use the PCMG construction of the coarser grids. 1477 1478 .seealso: PCMGGetGalerkin(), PCMGGalerkinType 1479 1480 @*/ 1481 PetscErrorCode PCMGSetGalerkin(PC pc,PCMGGalerkinType use) 1482 { 1483 PetscErrorCode ierr; 1484 1485 PetscFunctionBegin; 1486 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1487 ierr = PetscTryMethod(pc,"PCMGSetGalerkin_C",(PC,PCMGGalerkinType),(pc,use));CHKERRQ(ierr); 1488 PetscFunctionReturn(0); 1489 } 1490 1491 /*@ 1492 PCMGGetGalerkin - Checks if Galerkin multigrid is being used, i.e. 1493 A_i-1 = r_i * A_i * p_i 1494 1495 Not Collective 1496 1497 Input Parameter: 1498 . pc - the multigrid context 1499 1500 Output Parameter: 1501 . galerkin - one of PC_MG_GALERKIN_BOTH,PC_MG_GALERKIN_PMAT,PC_MG_GALERKIN_MAT, PC_MG_GALERKIN_NONE, or PC_MG_GALERKIN_EXTERNAL 1502 1503 Level: intermediate 1504 1505 .seealso: PCMGSetGalerkin(), PCMGGalerkinType 1506 1507 @*/ 1508 PetscErrorCode PCMGGetGalerkin(PC pc,PCMGGalerkinType *galerkin) 1509 { 1510 PC_MG *mg = (PC_MG*)pc->data; 1511 1512 PetscFunctionBegin; 1513 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1514 *galerkin = mg->galerkin; 1515 PetscFunctionReturn(0); 1516 } 1517 1518 PetscErrorCode PCMGSetAdaptInterpolation_MG(PC pc, PetscBool adapt) 1519 { 1520 PC_MG *mg = (PC_MG *) pc->data; 1521 1522 PetscFunctionBegin; 1523 mg->adaptInterpolation = adapt; 1524 PetscFunctionReturn(0); 1525 } 1526 1527 PetscErrorCode PCMGGetAdaptInterpolation_MG(PC pc, PetscBool *adapt) 1528 { 1529 PC_MG *mg = (PC_MG *) pc->data; 1530 1531 PetscFunctionBegin; 1532 *adapt = mg->adaptInterpolation; 1533 PetscFunctionReturn(0); 1534 } 1535 1536 PetscErrorCode PCMGSetAdaptCR_MG(PC pc, PetscBool cr) 1537 { 1538 PC_MG *mg = (PC_MG *) pc->data; 1539 1540 PetscFunctionBegin; 1541 mg->compatibleRelaxation = cr; 1542 PetscFunctionReturn(0); 1543 } 1544 1545 PetscErrorCode PCMGGetAdaptCR_MG(PC pc, PetscBool *cr) 1546 { 1547 PC_MG *mg = (PC_MG *) pc->data; 1548 1549 PetscFunctionBegin; 1550 *cr = mg->compatibleRelaxation; 1551 PetscFunctionReturn(0); 1552 } 1553 1554 /*@ 1555 PCMGSetAdaptInterpolation - Adapt the interpolator based upon a vector space which should be accurately captured by the next coarser mesh, and thus accurately interpolated. 1556 1557 Logically Collective on PC 1558 1559 Input Parameters: 1560 + pc - the multigrid context 1561 - adapt - flag for adaptation of the interpolator 1562 1563 Options Database Keys: 1564 + -pc_mg_adapt_interp - Turn on adaptation 1565 . -pc_mg_adapt_interp_n <int> - The number of modes to use, should be divisible by dimension 1566 - -pc_mg_adapt_interp_coarse_space <type> - The type of coarse space: polynomial, harmonic, eigenvector, generalized_eigenvector 1567 1568 Level: intermediate 1569 1570 .keywords: MG, set, Galerkin 1571 .seealso: PCMGGetAdaptInterpolation(), PCMGSetGalerkin() 1572 @*/ 1573 PetscErrorCode PCMGSetAdaptInterpolation(PC pc, PetscBool adapt) 1574 { 1575 PetscErrorCode ierr; 1576 1577 PetscFunctionBegin; 1578 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1579 ierr = PetscTryMethod(pc,"PCMGSetAdaptInterpolation_C",(PC,PetscBool),(pc,adapt));CHKERRQ(ierr); 1580 PetscFunctionReturn(0); 1581 } 1582 1583 /*@ 1584 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. 1585 1586 Not collective 1587 1588 Input Parameter: 1589 . pc - the multigrid context 1590 1591 Output Parameter: 1592 . adapt - flag for adaptation of the interpolator 1593 1594 Level: intermediate 1595 1596 .keywords: MG, set, Galerkin 1597 .seealso: PCMGSetAdaptInterpolation(), PCMGSetGalerkin() 1598 @*/ 1599 PetscErrorCode PCMGGetAdaptInterpolation(PC pc, PetscBool *adapt) 1600 { 1601 PetscErrorCode ierr; 1602 1603 PetscFunctionBegin; 1604 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1605 PetscValidPointer(adapt, 2); 1606 ierr = PetscUseMethod(pc,"PCMGGetAdaptInterpolation_C",(PC,PetscBool*),(pc,adapt));CHKERRQ(ierr); 1607 PetscFunctionReturn(0); 1608 } 1609 1610 /*@ 1611 PCMGSetAdaptCR - Monitor the coarse space quality using an auxiliary solve with compatible relaxation. 1612 1613 Logically Collective on PC 1614 1615 Input Parameters: 1616 + pc - the multigrid context 1617 - cr - flag for compatible relaxation 1618 1619 Options Database Keys: 1620 . -pc_mg_adapt_cr - Turn on compatible relaxation 1621 1622 Level: intermediate 1623 1624 .keywords: MG, set, Galerkin 1625 .seealso: PCMGGetAdaptCR(), PCMGSetAdaptInterpolation(), PCMGSetGalerkin() 1626 @*/ 1627 PetscErrorCode PCMGSetAdaptCR(PC pc, PetscBool cr) 1628 { 1629 PetscErrorCode ierr; 1630 1631 PetscFunctionBegin; 1632 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1633 ierr = PetscTryMethod(pc,"PCMGSetAdaptCR_C",(PC,PetscBool),(pc,cr));CHKERRQ(ierr); 1634 PetscFunctionReturn(0); 1635 } 1636 1637 /*@ 1638 PCMGGetAdaptCR - Get the flag to monitor coarse space quality using an auxiliary solve with compatible relaxation. 1639 1640 Not collective 1641 1642 Input Parameter: 1643 . pc - the multigrid context 1644 1645 Output Parameter: 1646 . cr - flag for compatible relaxaion 1647 1648 Level: intermediate 1649 1650 .keywords: MG, set, Galerkin 1651 .seealso: PCMGSetAdaptCR(), PCMGGetAdaptInterpolation(), PCMGSetGalerkin() 1652 @*/ 1653 PetscErrorCode PCMGGetAdaptCR(PC pc, PetscBool *cr) 1654 { 1655 PetscErrorCode ierr; 1656 1657 PetscFunctionBegin; 1658 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1659 PetscValidPointer(cr, 2); 1660 ierr = PetscUseMethod(pc,"PCMGGetAdaptCR_C",(PC,PetscBool*),(pc,cr));CHKERRQ(ierr); 1661 PetscFunctionReturn(0); 1662 } 1663 1664 /*@ 1665 PCMGSetNumberSmooth - Sets the number of pre and post-smoothing steps to use 1666 on all levels. Use PCMGDistinctSmoothUp() to create separate up and down smoothers if you want different numbers of 1667 pre- and post-smoothing steps. 1668 1669 Logically Collective on PC 1670 1671 Input Parameters: 1672 + mg - the multigrid context 1673 - n - the number of smoothing steps 1674 1675 Options Database Key: 1676 . -mg_levels_ksp_max_it <n> - Sets number of pre and post-smoothing steps 1677 1678 Level: advanced 1679 1680 Notes: 1681 this does not set a value on the coarsest grid, since we assume that 1682 there is no separate smooth up on the coarsest grid. 1683 1684 .seealso: PCMGSetDistinctSmoothUp() 1685 @*/ 1686 PetscErrorCode PCMGSetNumberSmooth(PC pc,PetscInt n) 1687 { 1688 PC_MG *mg = (PC_MG*)pc->data; 1689 PC_MG_Levels **mglevels = mg->levels; 1690 PetscErrorCode ierr; 1691 PetscInt i,levels; 1692 1693 PetscFunctionBegin; 1694 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1695 PetscValidLogicalCollectiveInt(pc,n,2); 1696 if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Must set MG levels with PCMGSetLevels() before calling"); 1697 levels = mglevels[0]->levels; 1698 1699 for (i=1; i<levels; i++) { 1700 ierr = KSPSetTolerances(mglevels[i]->smoothu,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,n);CHKERRQ(ierr); 1701 ierr = KSPSetTolerances(mglevels[i]->smoothd,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,n);CHKERRQ(ierr); 1702 mg->default_smoothu = n; 1703 mg->default_smoothd = n; 1704 } 1705 PetscFunctionReturn(0); 1706 } 1707 1708 /*@ 1709 PCMGSetDistinctSmoothUp - sets the up (post) smoother to be a separate KSP from the down (pre) smoother on all levels 1710 and adds the suffix _up to the options name 1711 1712 Logically Collective on PC 1713 1714 Input Parameters: 1715 . pc - the preconditioner context 1716 1717 Options Database Key: 1718 . -pc_mg_distinct_smoothup 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: PCMGSetNumberSmooth() 1727 @*/ 1728 PetscErrorCode PCMGSetDistinctSmoothUp(PC pc) 1729 { 1730 PC_MG *mg = (PC_MG*)pc->data; 1731 PC_MG_Levels **mglevels = mg->levels; 1732 PetscErrorCode ierr; 1733 PetscInt i,levels; 1734 KSP subksp; 1735 1736 PetscFunctionBegin; 1737 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1738 if (!mglevels) SETERRQ(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 const char *prefix = NULL; 1743 /* make sure smoother up and down are different */ 1744 ierr = PCMGGetSmootherUp(pc,i,&subksp);CHKERRQ(ierr); 1745 ierr = KSPGetOptionsPrefix(mglevels[i]->smoothd,&prefix);CHKERRQ(ierr); 1746 ierr = KSPSetOptionsPrefix(subksp,prefix);CHKERRQ(ierr); 1747 ierr = KSPAppendOptionsPrefix(subksp,"up_");CHKERRQ(ierr); 1748 } 1749 PetscFunctionReturn(0); 1750 } 1751 1752 /* No new matrices are created, and the coarse operator matrices are the references to the original ones */ 1753 PetscErrorCode PCGetInterpolations_MG(PC pc,PetscInt *num_levels,Mat *interpolations[]) 1754 { 1755 PC_MG *mg = (PC_MG*)pc->data; 1756 PC_MG_Levels **mglevels = mg->levels; 1757 Mat *mat; 1758 PetscInt l; 1759 PetscErrorCode ierr; 1760 1761 PetscFunctionBegin; 1762 if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling"); 1763 ierr = PetscMalloc1(mg->nlevels,&mat);CHKERRQ(ierr); 1764 for (l=1; l< mg->nlevels; l++) { 1765 mat[l-1] = mglevels[l]->interpolate; 1766 ierr = PetscObjectReference((PetscObject)mat[l-1]);CHKERRQ(ierr); 1767 } 1768 *num_levels = mg->nlevels; 1769 *interpolations = mat; 1770 PetscFunctionReturn(0); 1771 } 1772 1773 /* No new matrices are created, and the coarse operator matrices are the references to the original ones */ 1774 PetscErrorCode PCGetCoarseOperators_MG(PC pc,PetscInt *num_levels,Mat *coarseOperators[]) 1775 { 1776 PC_MG *mg = (PC_MG*)pc->data; 1777 PC_MG_Levels **mglevels = mg->levels; 1778 PetscInt l; 1779 Mat *mat; 1780 PetscErrorCode ierr; 1781 1782 PetscFunctionBegin; 1783 if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling"); 1784 ierr = PetscMalloc1(mg->nlevels,&mat);CHKERRQ(ierr); 1785 for (l=0; l<mg->nlevels-1; l++) { 1786 ierr = KSPGetOperators(mglevels[l]->smoothd,NULL,&(mat[l]));CHKERRQ(ierr); 1787 ierr = PetscObjectReference((PetscObject)mat[l]);CHKERRQ(ierr); 1788 } 1789 *num_levels = mg->nlevels; 1790 *coarseOperators = mat; 1791 PetscFunctionReturn(0); 1792 } 1793 1794 /*@C 1795 PCMGRegisterCoarseSpaceConstructor - Adds a method to the PCMG package for coarse space construction. 1796 1797 Not collective 1798 1799 Input Parameters: 1800 + name - name of the constructor 1801 - function - constructor routine 1802 1803 Notes: 1804 Calling sequence for the routine: 1805 $ my_csp(PC pc, PetscInt l, DM dm, KSP smooth, PetscInt Nc, const Vec initGuess[], Vec **coarseSp) 1806 $ pc - The PC object 1807 $ l - The multigrid level, 0 is the coarse level 1808 $ dm - The DM for this level 1809 $ smooth - The level smoother 1810 $ Nc - The size of the coarse space 1811 $ initGuess - Basis for an initial guess for the space 1812 $ coarseSp - A basis for the computed coarse space 1813 1814 Level: advanced 1815 1816 .seealso: PCMGGetCoarseSpaceConstructor(), PCRegister() 1817 @*/ 1818 PetscErrorCode PCMGRegisterCoarseSpaceConstructor(const char name[], PetscErrorCode (*function)(PC, PetscInt, DM, KSP, PetscInt, const Vec[], Vec **)) 1819 { 1820 PetscErrorCode ierr; 1821 1822 PetscFunctionBegin; 1823 ierr = PCInitializePackage();CHKERRQ(ierr); 1824 ierr = PetscFunctionListAdd(&PCMGCoarseList,name,function);CHKERRQ(ierr); 1825 PetscFunctionReturn(0); 1826 } 1827 1828 /*@C 1829 PCMGGetCoarseSpaceConstructor - Returns the given coarse space construction method. 1830 1831 Not collective 1832 1833 Input Parameter: 1834 . name - name of the constructor 1835 1836 Output Parameter: 1837 . function - constructor routine 1838 1839 Notes: 1840 Calling sequence for the routine: 1841 $ my_csp(PC pc, PetscInt l, DM dm, KSP smooth, PetscInt Nc, const Vec initGuess[], Vec **coarseSp) 1842 $ pc - The PC object 1843 $ l - The multigrid level, 0 is the coarse level 1844 $ dm - The DM for this level 1845 $ smooth - The level smoother 1846 $ Nc - The size of the coarse space 1847 $ initGuess - Basis for an initial guess for the space 1848 $ coarseSp - A basis for the computed coarse space 1849 1850 Level: advanced 1851 1852 .seealso: PCMGRegisterCoarseSpaceConstructor(), PCRegister() 1853 @*/ 1854 PetscErrorCode PCMGGetCoarseSpaceConstructor(const char name[], PetscErrorCode (**function)(PC, PetscInt, DM, KSP, PetscInt, const Vec[], Vec **)) 1855 { 1856 PetscErrorCode ierr; 1857 1858 PetscFunctionBegin; 1859 ierr = PetscFunctionListFind(PCMGCoarseList,name,function);CHKERRQ(ierr); 1860 PetscFunctionReturn(0); 1861 } 1862 1863 /* ----------------------------------------------------------------------------------------*/ 1864 1865 /*MC 1866 PCMG - Use multigrid preconditioning. This preconditioner requires you provide additional 1867 information about the coarser grid matrices and restriction/interpolation operators. 1868 1869 Options Database Keys: 1870 + -pc_mg_levels <nlevels> - number of levels including finest 1871 . -pc_mg_cycle_type <v,w> - provide the cycle desired 1872 . -pc_mg_type <additive,multiplicative,full,kaskade> - multiplicative is the default 1873 . -pc_mg_log - log information about time spent on each level of the solver 1874 . -pc_mg_distinct_smoothup - configure up (after interpolation) and down (before restriction) smoothers separately (with different options prefixes) 1875 . -pc_mg_galerkin <both,pmat,mat,none> - use Galerkin process to compute coarser operators, i.e. Acoarse = R A R' 1876 . -pc_mg_multiplicative_cycles - number of cycles to use as the preconditioner (defaults to 1) 1877 . -pc_mg_dump_matlab - dumps the matrices for each level and the restriction/interpolation matrices 1878 to the Socket viewer for reading from MATLAB. 1879 - -pc_mg_dump_binary - dumps the matrices for each level and the restriction/interpolation matrices 1880 to the binary output file called binaryoutput 1881 1882 Notes: 1883 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 1884 1885 When run with a single level the smoother options are used on that level NOT the coarse grid solver options 1886 1887 When run with KSPRICHARDSON the convergence test changes slightly if monitor is turned on. The iteration count may change slightly. This 1888 is because without monitoring the residual norm is computed WITHIN each multigrid cycle on the finest level after the pre-smoothing 1889 (because the residual has just been computed for the multigrid algorithm and is hence available for free) while with monitoring the 1890 residual is computed at the end of each cycle. 1891 1892 Level: intermediate 1893 1894 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCMGType, PCEXOTIC, PCGAMG, PCML, PCHYPRE 1895 PCMGSetLevels(), PCMGGetLevels(), PCMGSetType(), PCMGSetCycleType(), 1896 PCMGSetDistinctSmoothUp(), PCMGGetCoarseSolve(), PCMGSetResidual(), PCMGSetInterpolation(), 1897 PCMGSetRestriction(), PCMGGetSmoother(), PCMGGetSmootherUp(), PCMGGetSmootherDown(), 1898 PCMGSetCycleTypeOnLevel(), PCMGSetRhs(), PCMGSetX(), PCMGSetR() 1899 M*/ 1900 1901 PETSC_EXTERN PetscErrorCode PCCreate_MG(PC pc) 1902 { 1903 PC_MG *mg; 1904 PetscErrorCode ierr; 1905 1906 PetscFunctionBegin; 1907 ierr = PetscNewLog(pc,&mg);CHKERRQ(ierr); 1908 pc->data = (void*)mg; 1909 mg->nlevels = -1; 1910 mg->am = PC_MG_MULTIPLICATIVE; 1911 mg->galerkin = PC_MG_GALERKIN_NONE; 1912 mg->adaptInterpolation = PETSC_FALSE; 1913 mg->Nc = -1; 1914 mg->eigenvalue = -1; 1915 1916 pc->useAmat = PETSC_TRUE; 1917 1918 pc->ops->apply = PCApply_MG; 1919 pc->ops->applytranspose = PCApplyTranspose_MG; 1920 pc->ops->matapply = PCMatApply_MG; 1921 pc->ops->setup = PCSetUp_MG; 1922 pc->ops->reset = PCReset_MG; 1923 pc->ops->destroy = PCDestroy_MG; 1924 pc->ops->setfromoptions = PCSetFromOptions_MG; 1925 pc->ops->view = PCView_MG; 1926 1927 ierr = PetscObjectComposedDataRegister(&mg->eigenvalue);CHKERRQ(ierr); 1928 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCMGSetGalerkin_C",PCMGSetGalerkin_MG);CHKERRQ(ierr); 1929 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCMGGetLevels_C",PCMGGetLevels_MG);CHKERRQ(ierr); 1930 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCMGSetLevels_C",PCMGSetLevels_MG);CHKERRQ(ierr); 1931 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGetInterpolations_C",PCGetInterpolations_MG);CHKERRQ(ierr); 1932 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGetCoarseOperators_C",PCGetCoarseOperators_MG);CHKERRQ(ierr); 1933 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCMGSetAdaptInterpolation_C",PCMGSetAdaptInterpolation_MG);CHKERRQ(ierr); 1934 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCMGGetAdaptInterpolation_C",PCMGGetAdaptInterpolation_MG);CHKERRQ(ierr); 1935 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCMGSetAdaptCR_C",PCMGSetAdaptCR_MG);CHKERRQ(ierr); 1936 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCMGGetAdaptCR_C",PCMGGetAdaptCR_MG);CHKERRQ(ierr); 1937 PetscFunctionReturn(0); 1938 } 1939